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Preface 


The Algorithms such as SVD, Eigen decomposition, Gaussian Mixture 
Model, PSO, Ant Colony etc. are scattered in different fields. There is the 
need to collect all such algorithms for guick reference. Also there is the need 
to view such algorithms in application point of view. This Book attempts to 
satisfy the above reguirement. Also the algorithms are made clear using 
MATLAB programs. This book will be useful for the Beginners Research 
scholars and Students who are doing research work on practical applications 
of Digital Signal Processing using MATLAB. 
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Chapter 1 


ARTIFICIAL INTELLIGENCE 
Algorithm Collections 


1. PARTICLE SWARM ALGORITHM 


Consider the two swarms flying in the sky, trying to reach the particular 
destination. Swarms based on their individual experience choose the proper 
path to reach the particular destination. Apart from their individual 
decisions, decisions about the optimal path are taken based on their 
neighbor’s decision and hence they are able to reach their destination faster. 
The mathematical model for the above mentioned behavior of the swarm is 
being used in the optimization technique as the Particle Swarm Optimization 
Algorithm (PSO). 

For example, let us consider the two variables ‘x’ and ‘y’ as the two 
swarms. They are flying in the sky to reach the particular destination (1.е.) 
they continuously change their values to minimize the function (х-10)-Қу- 
5)’. Final value for ‘x’ and “у” are 10.1165 and 5 respectively after 100 
iterations. 

The Figure 1-1 gives the closed look of how the values of x and y are 
changing along with the function value to be minimized. The minimization 
function value reached almost zero within 35 iterations. Figure 1-2 shows 
the zoomed version to show how the position of x and y are varying until 
they reach the steady state. 
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Figure 1-1. PSO Example zoomed version 
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Figure 1-2. PSO Example 


1.1 How are the Values of ‘x and y’ are Updated 
in Every Iteration? 


The vector representation for updating the values for x and y is given in 
Figure 1-3. Let the position of the swarms be at ‘a’ and ‘b’ respectively as 
shown in the figure. Both are trying to reach the position ‘e’. Let ‘a’ decides 
to move towards “с” and ‘b’ decides to move towards “а”. 

The distance between the position “с” and “е” is greater than the distance 
between 'd' and ‘e’. so based on the neighbor's decision position ‘d’ is 
treated as the common position decided by both ‘a’ and ‘b’. (ie) the position 
“с? is the individual decision taken by ‘a’, position ‘d’ is the individual 
decision taken by ‘b’ and the position “а” is the common position decided by 
both ‘a’ and ‘b’. 
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Figure 1-3. Vector Representation of PSO Algorithm 


‘a’ based on the above knowledge, finally decides to move towards the 
position ‘g’ as the linear combination of ‘oa’ , ‘ac’ and ‘ad’. [As ‘d’ is the 
common position decided].The linear combination of ‘oa’ and scaled ‘ac’ 
(ie) “аР is the vector ‘of’. The vector ‘of? combined with vector ‘fg’ (ie) 
scaled version of ‘ad’ to get ‘og’ and hence final position decided by ‘a’ is 


Le 


р’. 
Similarly, ‘b’ decides the position ‘h’ as the final position. It is the linear 
combination of ‘ob’ and ‘bh’(ie) scaled version of ‘bd’. Note as ‘d’ is the 
common position decided by ‘a’ and ‘b’, the final position is decided by 
linear combinations of two vectors alone. 

Thus finally the swarms ‘a’ and ‘b’ moves towards the position ‘g’ and 
‘h’ respectively for reaching the final destination position ‘e’. The swarm ‘a’ 
and ‘b’ randomly select scaling value for linear combination. Note that ‘oa’ 
and ‘ob’ are scaled with 1 (ie) actual values are used without scaling. Thus 
the decision of the swarm ‘a’ to reach “е” is decided by its own intuition 
along with its neighbor’s intuition. 

Now let us consider three swarms (A,B,C) are trying to reach the 
particular destination point *D'. A decides A’, B decides В’ and С decides 
С” as the next position. Let the distance between the B' and D is less 
compared with A’D and С” and hence, В” is treated as the global decision 
point to reach the destination faster. 

Thus the final decision taken by A is to move to the point, which is the 
linear combination of OA, AA' and AB'. Similarly the final decision taken 


4 Chapter 1 


by Bis to move the point which is the linear combination of OB, BB”. The 
final decision taken by C is to move the point which is the linear 
combination of ОС, СС? and CB”. 


1.2 PSO Algorithm to Maximize the Function F (X, Y, Z) 


1. Initialize the values for initial position a, b, c, d, e 

2. Initialize the next positions decided by the individual swarms as a’, b’, c 
d' and e 

3. Global decision regarding the next position is computed as follows. 
Compute f (a’, b, с, d, е), f (a, b’, c, d, e), f(a, b, с”, а, е), f (a, b, c, d^, е) 
and f (a, b, c, d, е”). Find minimum among the computed values. If f (a’, 
b, c, d, e) is minimum among all, the global position decided regarding 
the next position is a’. Similarly If f (a, b’, c, d, e) is minimum among all, 
b' is decided as the global position regarding the next position to be 
shifted and so on. Let the selected global position is represented ad 
‘global’ 

4. Next value for a is computed as the linear combination of ‘a’ , (a’-a) and 
(global-a) (ie) 


, 


nexta = a+ СІ * RAND * (a’ —a) + C2 * RAND * (global —a ) 
nextb = b+ СІ * RAND * (b' —b) + C2 * RAND * (global —b) 
nextc = c+ СІ * RAND * (c? —c) + C2 * RAND * (global —c ) 
nextd = d+ СІ * RAND * (d' ~) + C2 * RAND * (global —d ) 
nexte = e+ СІ * RAND * (e’ —e) + C2 * RAND * (global —e ) 














5. Change the current value for a, b, c, d and e as nexta, nextb, nextc, nextd 
and nexte 

6. If f (nexta, b, c, d, e) is less than f (a’, b, c, d, e) then update the value for 
а” as nexta, otherwise a’ is not changed. 


If f (a, nextb, c, d, e) is less than f (a, b’, c, d, e) then update the value for 
b’ as nextb, otherwise b’ is not changed 
If f (a, b, nextc, d, e) is less than f (a, b, с”, d, e) then update the value 
for с” as nextc, otherwise с” is not changed 
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If f (a, b, с, nextd, e) is less than f (a, b, с, 4”, е) then update the value 
for d’ as nextd, otherwise d’ is not changed 


If f (a, b, c, d, nexte) is less than f(a, b, c, d, е”) then update the value 
for e' as nexte, otherwise e' is not changed 


7. Repeat the steps 3 to 6 for much iteration to reach the final decision. 


The values for ‘cl’,’c2’ are decided based on the weightage given to 
individual decision and global decision respectively. 

Let Aa(t) is the change in the value for updating the value for ‘a’ in tth 
iteration, then nexta at (t-1)th iteration can be computed using the following 
formula. This is considered as the velocity for updating the position of the 
swarm in every iteration. 


nexta (t+1) = a (t) + Aa(t+1) 

where 

Aa(t+1)= cl * rand * (a*—a ) + c2 * rand * ( global —a) + 
w(t)* Aa(t) 





‘w (t) is the weight at t” iteration. The value for ‘w’ is adjusted at every 
iteration as given below, where ‘iter’ is total number of iteration used. 


w(t+1)=w(t)-t*w(t)/(iter). 


Decision taken in the previous iteration is also used for deciding the next 
position to be shifted by the swarm. But as iteration increases, the 
contribution of the previous decision is decreases and finally reaches zero in 
the final iteration. 
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psogv .m 


function [value]-psogv(fun,range,ITER) 
Yopsogv.m 

%Particle swarm algorithm for maximizing the function fun with two variables x 
%and y. 

%Syntax 
%[value]=psogv(fun,range, ITER) 
%example 

%fun="f1' 

%create the function fun.m 

%function [res ]=fun(x,y) 
Yres=sin(x)+cos(x); 

%range=[-pi pi;-pi рі]; 

%ITER is the total number of Iteration 
error=[]; 

vell=[]; 

vel2=[]; 


%Intialize the swarm position 
swarm=[]; 
x(1)=rand*range(1,2)+range(1, 1); 
y(1)=rand*range(2,2)+range(2, 1); 
x(2)=rand*range(1,2)+range(1, 1); 
y(2)=rand*range(2,2)+range(2, 1); 





“%Intialize weight 

w-l; 

cl=2; 

c2=2; 

%Initialize the velocity 

v1=0;%velocity for x 

v2=0;%velocity for y 

for i=1:1:ITER 

[p.g]=min([f1 (fun,x(2),y(1)) Р (fun,x(1),y(2)))); 


if (q==1) 
capture=x(2); 
else 


capture=y(2); 
end 
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Continued... 


vl=w*v1+c1 *rand*(x(2)-x(1))4 
v2=w*v2+c1 *rand*(y(2)-y(1))- 
vell=[vel1 v1]; 
vel2=[vel2 v2]; 








Youpdating x(1) and y(1) 
х(1)=х(1)+у1; 
y()=y(1)+v2; 





%updating x(2) and y(2) 


-c2*rand*(capture-x(1)); 
-c2*rand*(capture-y(1)); 


ICI (fun,x(2),y(1)))<=(F1 fun;x(1),y(1)))) 


х(2)=х(2); 
else 
x(2)=x(1); 


end; 


II (оп,х(1),у(2)))<=(# (fun,x(1),y(1))) 


у(2)=у(2); 

else 

y(2)=y(1); 

end 

error=[error f1(fun,x(2),y(2))]; 
w-w-w*i/ITER; 
swarm=[swarm;x(2) y(2)]; 
subplot(3,1,3) 
plot(error,'-') 
title('Error(vs) Iteration’); 
subplot(3, 1,1) 
plot(swarm(:,1),'-" 

title(x (vs) Iteration'); 
subplot(3,1,2) 
plot(swarm(:,2),'-') 

title(y (vs) Iteration’); 
pause(0.2) 

end 

value=[x(2);y(2)]; 





fl.m 


function [res]=f1(fun,x,y); 
s=streat(fun,'(x,y)'); 
res=eval(s); 
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1.4 Program Illustration 


Following the sample results obtained after the execution of the program 
psogv.m for maximizing the function 'f1.m" 
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2. GENETIC ALGORITHM 


A basic element of the Biological Genetics is the chromosomes. 
Chromosomes cross over each other. Mutate itself and new set of 
chromosomes is generated. Based on the reguirement, some of the 
chromosomes survive. This is the cycle of one generation in Biological 
Genetics. The above process is repeated for many generations and finally 
best set of chromosomes based on the reguirement will be available. This is 
the natural process of Biological Genetics. The Mathematical algorithm 
eguivalent to the above behavior used as the optimization technigue is called 
as Artificial Genetic Algorithm. 

Let us consider the problem for maximizing the function f(x) subject to 
the constraint x varies from ‘m’ to ‘n’. The function f(x) is called fitness 
function. Initial population of chromosomes is generated randomly. (1.e.) the 
values for the variable ‘x’ are selected randomly between the range ‘m’ to 
‘n’. Let the values be хі, Хә.....хі, where ‘L’ is the population size. Note that 
they are called as chromosomes in Biological context. 

The Genetic operations like Cross over and Mutation are performed to 
obtain ‘2*L’ chromosomes as described below. 

Two chromosomes of the current population is randomly selected (ie) 
select two numbers from the current population. Cross over operation 
generates another two numbers у1 and y2 using the selected numbers. Let 
the randomly selected numbers be x3 and x9. Yı is computed as r*x3+(1- 
r)*x9. Similarly ух is computed as (1-r)*x3*r*xo, where ‘r’ is the random 
number generated between 0 tol. 

The same operation is repeated ‘L’ times to get '2*L" newly generated 
chromosomes. Mutation operation is performed for the obtained 
chromosomes to generate ‘2*L’ mutated chromosomes. For instance the 
generated number ‘yı? is mutated to give zı mathematically computed as 
r*y, where гү is the random number generated. Thus the new set of 
chromosomes after crossover and Mutation are obtained as [zi 72 Z3 ...Z21]. 

Among the ‘2L’ values generated after genetic operations, ‘L’ values аге 
selected based on Roulette Wheel selection. 
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2.1 Roulette Wheel Selection Rule 


Consider the wheel partitioned with different sectors as shown in the 
Figure 1-4. Let the pointer *p' be in the fixed position and wheel is pivoted 
such that the wheel can be rotated freely. This is the Roulette wheel setup. 
Wheel is rotated and allowed to settle down. 

The sector pointed by the pointer after settling is selected. Thus the 
selection of the particular sector among the available sectors are done using 
Roulette wheel selection rule. 

In Genetic flow ‘L’ values from ‘2L’ values obtained after cross over and 
mutation are selected by simulating the roulette wheel mathematically. 
Roulette wheel is formed with ‘2L’ sectors with area of each sector is 
proportional to f(z1), f(z;) f(z3)... and {22 ) respectively, where ‘f(x)’ is the 
fitness function as described above. They are arranged in the row to form the 
fitness vector as [f(zi), f(Z2) f(z3)...f(zoi )]. Fitness vector is normalized to 
form Normalized fitness vector as [f,(z1), fi(Z2) fr(Z3)...fn(Zar )], so that sum 
of the Normalized fitness values become 1 (1.е.) normalized fitness value of 
f(z,) is computed as (21) = 1) / [ 21) + Қо)- Қаз)- f(z4)... KZ21)]. 
Similarly normalized fitness value is computed for others. 

Cumulative distribution of the obtained normalized fitness vector is 
obtained as 


[f (zi) f(zi)* (22) fi(zi)* f(z5)* (23)... 1]. 


Generating the random number 'r' simulates rotation of the Roulette 
Wheel. 

Compare the generated random number with the elements of the 
cumulative distribution vector. If ‘r< f.(z ) and ‘r > 0°, the number ‘z’ is 
selected for the next generation. Similarly if ‘r< f,(z,)+ f,(Z2)’ and ‘r > f,(z1)’ 
the number ‘z?’ is selected for the next generation and so on. 


Figure 1-4. Roulette Wheel 
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The above operation defined as the simulation for rotating the roulette 
wheel and selecting the sector is repeated ‘L’ times for selecting ‘L’ values 
for the next generation. Note that the number corresponding to the big sector 
is having more chance for selection. 

The process of generating new set of numbers (ie) next generation 
population from existing set of numbers (ie) current generation population is 
repeated for many iteration. 

The Best value is chosen from the last generation corresponding to the 
maximum fitness value. 


2.2 Example 


Let us consider the optimization problem for maximizing the function f(x) = 
x+10*sin (5*x) +7*cos (4*x) +sin(x) ,subject to the constraint x varies from 
0 to 9. 

The numbers between 0 and 9 with the resolution of 0.001 are treated as 
the chromosomes used in the Genetic Algorithm. (ie) Float chromosomes are 
used in this example. Population size of 10 chromosomes is survived in 
every generation. Arithmetic cross over is used as the genetic operator. 
Mutation operation is not used in this example Roulette Wheel selection is 
made at every generation. Algorithm flow is terminated after attaining the 
maximum number of iterations. In this example Maximum number of 
iterations used is 100. 

The Best solution for the above problem is obtained in the thirteenth 
generation using Genetic algorithm as 4.165 and the corresponding 
fitness function f(x) is computed as 8.443. Note that Genetic algorithm 
ends up with local maxima as shown in the figure 1-5. This is the 
drawback of the Genetic Algorithm. When you run the algorithm again, it 
may end up with Global maxima. The chromosomes generated randomly 
during the first generation affects the best solution obtained using genetic 
algorithm. 

Best chromosome at every generation is collected. Best among the 
collection is treated as the final Best solution which maximizes the 
function f(x). 


2.2.1 M-program for genetic algorithm 


The Matlab program for obtaining the best solution for maximizing the 
function f(x) = x+10*sin (5*x) +7*cos (4*x) +sin(x) using Genetic 
Algorithm is given below. 
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geneticgv.m 


clear all, close all 
pop=0:0.001:9; 
pos=round(rand(1,10)*9000)+1; 
pop 1=pop(pos); 
BEST=[]; 
for iter=1:1:100 
coll=[]; 
col2=[]; 
for do=1:1:10 
rl=round(rand(1)*9)+1; 
r2=round(rand(1)*9)+1 ; 
r3=rand; 
vl1=r3*pop1(11)+(1-13)*pop1(12); 
v2=r3*pop1(12)+(1-13)*pop1(11); 
coll=[coll v1 v2]; 
end 
sect=fcn(col1)+abs(min(fcn(col1))); 
sect=sect/sum(sect); 
[u,v]J=min(sect); 
c=cumsum(sect); 
for i=1:1:10 
r=rand; 
cl=c-r; 
[p,g]=find(e1>=0); 
if(length(g)-=0) 
col2=[col2 col1(q(1))]; 
else 
col2=[col2 coll(v)]; 
end 
end 
pop 1=col2; 
s=fcn(pop); 
plot(pop,s) 
[u,v]=max(fen(pop1)); 
BEST=[BEST;pop1(v(1)) fcn(pop1(v(1)))]; 
hold on 
plot(pop1,fcn(pop1),T."); 
M(iter)=getframe; 
pause(0.3) 
hold off 
[iter pop1(v(1)) fen(popl(v(1))] 
end 
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рор1=со12; 
s=fcn(pop); 
plot(pop,s) 
[u,v]J=max(fcn(pop1)); 
BEST=[BEST;pop1(v(1)) fcn(pop1(v(1)))]; 
hold on 
plot(pop1,fcn(pop1),T."); 
M(iter)=getframe; 
pause(0.3) 
hold off 
[iter pop 1(v(1)) fen(popl(v(1))] 
end 
for =1:1:14 
D(:,:,:,1)=M(i).cdata; 
end 
figure 
imshow(M(1).cdata) 
figure 
imshow(M(4).cdata) 
figure 
imshow(M(10).cdata) 
figure 
imshow(M(30).cdata) 





fcn.m 


function [res]=fcn(x) 
res=x+10*sin(5*x)+7*cos(4*x)+sin(x); 





Note that the m-file fen.m may be edited for changing the fitness function 


2.2.2 Program illustration 


The following is the sample results obtained during the execution of the 


program geneticgv.m 
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Figure 1-5. Convergence of population 


The graph is plotted between the variable ‘x’ and fitness function “(х)”. The 
values for the variable ‘x’ ranges from 0 to 9 with the resolution of 0.001. The 
dots on the graph are indicating the fitness value of the corresponding 
population in the particular generation. Note that in the first generation the 
populations are randomly distributed. In the thirteenth generation populations 
are converged to the particular value called as the best value, which is egual to 
4.165, and the corresponding fitness value f (x) is 8.443. 
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Classification of Genetic Operators 


In genetic algorithm, chromosomes are represented as follows 


> Float form 
The raw data itself acts as the chromosomes and can be used without 
without any modification. 


> Binary form 

Binary representation of the number is used as the chromosomes. Num- 
ber of Bits used to represent the raw data depends upon the resolu- 
tion required which is measure of accuracy in representation. 
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> Order based form 


Suppose the reguirement 15 to find the order in which sales man is 
traveling among the ‘n’ places such that total distance traveled by that 
person is minimum. In this case, the data is the order in which the characters 
are arranged. Example data = [a b e f g c]. This type of representation is 
called Order based form of representation. 

The tree diagram given above summarizes the different type of Genetic 
operators. The Examples for every operators are listed below 


2.3.1 Simple crossover 


Suppose the requirement is to maximize the function f (x, у, z) subject to the 
constraints x varies from 0.2 to 0.4, y varies from 0.3 to 0.6 and z varies 
from 0.1 to 0.4. 

Let the two chromosomes of the current population be Cl and C2 
respectively. 


Before crossover 
C1 = [0.31 0.45 0.11] and C2 = [0.25 0.32 0.3] (say). 
After crossover 


СІ = [0.3100 0.3200 0.3000] 
C2 = [0.2500 0.4500 0.1100] 


Simple cross over can also be applied for Binary representation as 
mentioned below 


Before Crossover 


C1 = [10101°, “01010°,*01011°] 
C2 = [11111°, 010117, “101007 


After Crossover 


C1 =[‘10101’,’01010’,’011007] 
C2 = [:11111,01011', 100117] 


Note that the crossover point is selected by combining all the binary string 
into the single string. 


2.3.2 Heuristic crossover 


Suppose the requirement is to maximize the function f (x, у, z) subject to the 
constraints x varies from 0.2 to 0.4, y varies from 0.3 to 0.6 and z varies 
from 0.1 to 0.4. 
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Let the two chromosomes of the current population be C1 and C2 
respectively. 


Before crossover 


СІ = [0.31 0.45 0.11] and C2 = [0.25 0.32 0.3] (say) and the 
corresponding fitness value (ie) f(x,y,z)=0.9 and 0.5 respectively 

Heuristic crossover identifies the best and worst chromosomes as 
mentioned below 

Best chromosome is one for which the fitness value is maximum. Worst 
chromosome is one for which the fitness value is minimum. In this example 
Best chromosome [Bt] is [0.31 0.45 0.11] and Worst Chromosome [Wt] is 
[0.25 0.32 0.3]. 


Bt = [0.31 0.45 0.11] 
Wt = [0.25 0.32 0.3] 


Bt chromosome 1s returned without disturbance (ie) After cross over, 
CI-Bt- [0.31 0.45 0.11] 


Second Chromosome after crossover is obtained as follows. 
1. Generate the random number ‘r’ 
C2=(Bt-Wt)*r+Bt 


2. Check whether C2 is within the required boundary values (i.e.) first value 
of the vector C2 is within the range from 0.2 to 0.4, second value is within 
the range from 0.3 to 0.6 and third value is within the range from 0.1 to 0.4. 
3. If the condition is satisfied, computed C2 is returned as the second 
chromosome after crossover. 

4. If the condition is not satisfied, repeat the steps 1 to 3. 

5. Finite attempts are made to obtain the modified chromosome after cross 
over as described above. (i.e) If the condition is not satisfied for ‘N’ attempts, 
C2=Wt is returned as the second chromosome after crossover. 


2.3.3 Arith crossover 


Let the two chromosomes of the current population be СІ and C2 
respectively. 
Before crossover 


C1 = [0.31 0.45 0.11] and C2 = [0.25 0.32 0.3] (say) 
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After Crossover 


Generate the random number r=0.3(say) 


Cl=r*C1+(1-r)*C2= [0.2680 0.3590 0.2430] 
C2=(1-r)*C1+r*C2= [0.2920 0.4110 0.1670 ] 


3. SIMULATED ANNEALING 


The process of heating the solid body to the high temperature and allowed to 
cool slowly is called Annealing. Annealing makes the particles of the solid 
material to reach the minimum energy state. This is due to the fact that when 
the solid body is heated to the very high temperature, the particles of the 
solid body are allowed to move freely and when it is cooled slowly, the 
particles are able to arrange itself so that the energy of the particles are made 
minimum. The mathematical equivalent of the thermodynamic annealing as 
described above is called simulated annealing. 

The energy of the particle in thermodynamic annealing process can be 
compared with the cost function to be minimized in optimization problem. 
The particles of the solid can be compared with the independent variables 
used in the minimization function. 

Initially the values assigned to the variables are randomly selected from 
the wide range of values. The cost function corresponding to the selected 
values are treated as the energy of the current state. Searching the values 
from the wide range of the values can be compared with the particles 
flowing in the solid body when it is kept in high temperature. 

The next energy state of the particles is obtained when the solid body is 
slowly cooled. This is equivalent to randomly selecting next set of the 
values. 

When the solid body is slowly cooled, the particles of the body try to 
reach the lower energy state. But as the temperature is high, random flow of 
the particles still continuous and hence there may be chance for the particles 
to reach higher energy state during this transition. Probability of reaching the 
higher energy state is inversely proportional to the temperature of the solid 
body at that instant. 

In the same fashion the values are randomly selected so that cost function 
of the currently selected random values is minimum compared with the 
previous cost function value. At the same time, the values corresponding to 
the higher cost function compared with the previous cost function are also 
selected with some probability. The probability depends upon the current 
simulated temperature ‘T’. If the temperature is large, probability of 
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selecting the values corresponding to higher energy levels are more. This 
process of selecting the values randomly is repeated for the finite number of 
iteration. The values obtained after the finite number of iteration can be 
assumed as the values with lowest energy state (i.e) lowest cost function 
Thus the simulated Annealing algorithm is summarized as follow. 


3.1 Simulated Annealing Algorithm 


The optimization problem is to estimate the best values for ‘x’, ‘y’ and 
‘z’ such that the cost function f (x, у, 2) is minimized. ‘x’ varies from 
*Xmin ‘tO ‘Xmax’. ‘У’ varies from ‘y’ varies “Ymin' to ‘Ymax .'7' varies from 
€ * € % 
Zmin to Zmax 


Step 1: Intialize the value the temperature ‘T’. 

Step 2: Randomly select the current values for the variables ‘x’, °y’ and ‘z’ 
from the range as defined above. Let it be ‘xe, ‘Ye and ‘zè 
respectively. 

Step 3: Compute the corresponding cost function value f (х, у, Zo). 

Step 4: Randomly select the next set of values for the variables ‘x’, °y’ and 
‘z’ from the range as defined above. Let it be ‘xy’, ‘Yn’ and ‘Zr 
respectively. 

Step 5: Compute the corresponding cost function value f (Xn, Yn, Zn). 

Step 6: If f (x, Yn, Zn)<= f (Xe, Ve, Zc), then the current values for the random 
variables хұ- Xn, Yc,- Yn and Z.=Zp. 

Step 7: If f(xy, Yn, Zn)>f(Xc, Ус, Zc), then the current values for the random 
variables хұ- Xn, Ус- Yn and 7-27, are assigned only when 


exp [(f (xo уг, Zo)- f(Xn, Yn, Zn))/T] > rand 


Note that when the temperature “Т” is less, the probability of selecting the 
new values as the current values is less. 


Step 8: Reduce the temperature Т=г*Т. r' is scaling factor varies from 
0 to 1. 

Step 9: Repeat STEP 3 to STEPS for n times until “Т” reduces to the 
particular percentage of initial value assigned to *T* 


3.2 Example 


Consider the optimization problem of estimating the best values for *x' such 
that the cost function f(x)» x+10*sin(5*x)+7*cos(4*x)+sin(x) is minimized 
and x varies from 0 to 5. 
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Figure 1-6 Illustration of simulated Annealing 


The figure 1-6 shows all possible cost function values corresponding to 
the ‘x’ values ranges from 0 to 5. The goal is to find the best value of x 
corresponding to global minima of the graph given above. Note that the 
global minima is approximately at x=0.8. 

Initially the random value is selected within the range (0 to 5). Let the 
selected value be 2.5 and the corresponding cost function is -3.4382. The 
variable ‘curval’ is assigned the value 2.5 and the variable ‘curcost’ is 
assigned the value -3.4382. Intialize the simulated temperature as T=1000. 

Let the next randomly selected value be 4.9 and the corresponding cost 
function is 3.2729. The variable ‘newval’ is assigned the value 4.9 and the 
*newcost' is assigned the value 3.2729. 

Note that ‘newcost’ value is higher than the “curcost* value. As 
‘newcost’> ‘curcost’, ‘newval’ is inferior compared with ‘curval’ in 
minimizing the cost function. As the temperature (T) is large and 
exp((curcost-newcost)/T)>rand, ‘newcost’ value is assigned to ‘curcost’ and 
‘newval’ is assigned to ‘curval’. Now the temperature “Т? is reduced by the 
factor 0.2171=1/log(100), where 100 is the number of iterations used. This is 
the thumb rule used in this example. This process is said to one complete 
iteration. Next randomly selected value is selected and the above described 
process is repeated for 100 iterations. 

The order in which the values are selected in the first 6 iterations is not 
moving towards the local minima point which can be noted from the figure 
1-6. This is due to the fact that the initial simulated temperature of the 
annealing process is high. 
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Randomly selected values for the varible 'x' and f(x) 
for the first 6 Iterations 
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Figurel-7 Illustration of Simulated Annealing 1 


As iteration increases the simulated temperature is decreased and the 
value selected for the variable ‘x is moving towards the global minima point 
as shown in the figure 1-7. 

Thus values assigned to the variable ‘x’ for the first few iterations 
(about 10 iterations in this example) is not really in decreasing order of 
f(x). This helps to search the complete range and hence helps in 
overcoming the problem of local minima. As iteration goes on increasing 
the values selected for the variable ‘x’ is moving towards the global 
minima which can be noted from the figure 1-8. 
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The value assianed to the variable 'x' and f(x) for the 1.10.20.30 Iterations 
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Figure 1-8. Illustration of Simulated Annealing 2 
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Figurel-9. Illustration of Simulated Annealing 3 
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3.3 M-program for Simulated Annealing 


Matlab program for minimizing the function f(x)= x+10*sin(5*x)+ 
7*cos(4*x)+sin(x),where x varies from 0 to 5. 





sagv.m 


x=0:1/300:5 
plot(x, minfcn(x)) 
hold 
curval=rand*5; 
curcost=minfcn(curval); 
T=1000; 
for 1=1:1:100 
newval=rand*5; 
newcost=minfcn(newval); 
if((newcost-curcost)<=0) 
curval=newval; 
curcost=minfcn(curval); 
else 
if(exp((curcost-newcost)/T)>rand) 
curval=newval; 
curcost=minfcn(curval); 
end 
end 
T=T/log(100) 
plot(curval,minfcn(curval),'r*"); 
CURRENTVAL {i}=curval; 
М {i}=minfcn(curval); 
pause(0.2) 
end 





minfcn.m 


function [res]=minfcn(x) 
res=x+10*sin(5*x)+7*cos(4*x)+sin(x) 
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4. BACK PROPAGATION NEURAL NETWORK 


The mathematical model of the Biological Neural Network is defined as 
Artificial Neural Network. One of the Neural Network models which are 
used almost in all the fields is Back propagation Neural Network. 

The model consists of layered architecture as shown in the figure 1-10. 

‘T is the Input layer “О” is the output layer and ‘H’ is the hidden layer. 
Every neuron of the input layer is connected to every neuron in the hidden 
layer. Similarly every neuron in the Hidden layer is connected with every 
neuron in the output layer. The connection is called weights. 

Number of neurons in the input layer is equal to the size of the input 
vector of the Neural Network. Similarly number of neurons in the output 
layer is equal to the size of the output vector of the Neural Network. Size of 
the hidden layer is optional and altered depends upon the requirement 

For every input vector, the corresponding output vector is computed as 
follows 


[hidden vector] = func ([Input vector]*[Weight Matrix 1] + [Bias1 |) 
[output vector] = func ([hidden vector]*[Weight Matrix 2 | + [Bias2 |) 


If the term [Input vector]*[Weight Matrix 1] value becomes zero, the 
output vector also becomes zero. This can be avoided by using Bias vectors. 
Similarly to make the value of the term “[Input vector]*[Weight Matrix 1] + 
[В1а$1]” always bounded, the function ‘func’ is used as mentioned in the 
equation given above. 


| O 


Figure 1-10. BPN Structure 
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commonly Usually used ‘func’ are 


logsig (x) = 1 


[1+exp(-x)] 


tansig(x) = 2 


[1+exp(-2x)] 


The reguirement is to find the common Weight Matrix and common Bias 
vector such that for the particular input vector, computed output vector must 
match with the expected target vector. The process of obtaining the Weight 
matrix and Bias vector is called training. 

Consider the architecture shown in the figure 1-9. Number of neurons in 
the input layer is 3 and number of neurons in the output layer is 2. Number of 
neurons in the hidden layer is 1.Тһе weight connecting the i^ neuron in the 
first layer and j neuron in the hidden layer is represented as Wj, The weight 
connecting i^ neuron in the hidden layer and j" neuron in the output layer is 
represented as Wj; . 

Let the input vector is represented as [i; 12 13]. Hidden layer output is 
represented as [h;] and output vector is represented as [o; 02]. The bias 
vector in the hidden layer is given as [b,] and the bias vector in the output 
layer is given as |р, b2]. Desired ouput vector is represented is [t1 t2] 


The vectors are related as follows. 
hl=funcl (Wifi | үу; 12 t W31*l3 t by ) 
01= func2 (wi *h; tbi) 

02- func2 (wi? *hi +b2) 

tl~=ol 

t2~=02 








4.1 Single Neuron Architecture 


Consider the single neuron input layer and single neuron output layer. 

Output= func (input * W +B) 

Desired ouput vector be represented as ‘Target’ 

Requirement is to find the optimal value of W so that Cost function = 
(Target-Output)? is reduced. The graph plotted between the weight vector 
‘W’ and the cost function is given in the figure 1-11. 

Optimum value of the weight vector is the vector corresponding to the 
point 1 which is global minima point, where the cost value is lowest. 
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cost function 


weight vector 


Figure 1-11 ANN Algorithm Illustration 


The slope computed at point 2 in the figure is negative. Suppose if the 
weight vector corresponding to the point 2 is treated as current weight vector 
assumed, the next best value (i.e) global point is at right side of it. Similarly 
if the slope is positive, next best value is left side of the current value. This 
can be seen from the graph. The slope at point 3 is a positive. The best value 
(i.e) global minimum is left of the current value. 

Thus to obtain the best value of the weight vector. Initialize the weight 
vector. Change the weight vector iteratively. Best weight vector is obtained 
after finite number of iteration. The change in the weight vector in every 
iteration is represented as ‘AW’. 





i.e.) W (0+1) =W (n) + AW (n+1) 


W(n+1) is the weight vector at (n+1)" iteration W(n) is the weight vector 
at n" iteration and AW(n+1) is the change in the weight vector at (п+1)" 
iteration. 

The sign and magnitude of the ‘AW’ depends upon the direction and 
magnitude of the slope of the cost function computed at the point 
corresponding to the current weight vector. 
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Thus the weight vector is adjusted as following 
New weight vector = Old weight vector - u*slope at the current cost value. 
Cost function(C) = (Target-Output)* 
=> C= (Target- func (input * W +B)) 


Slope is computed by differentiating the variable C with respect to W and 
is approximated as follows. 


2((Target- func (input * W +B))* (-input) 
=2 *error*input 


Thus W(n+1) =W(n) - 2*p*error*(-input) 





=> W(n+1) = W(n) + ү error * input 


This is back propagation algorithm because error is back propagated for 
adjusting the weights. The Back propagation algorithm for the sample 
architecture mentioned in figure 1-9 is given below. 


4.2 Algorithm 


Step 1: Initialize the Weight matrices with random numbers. 

Step 2: Initialize the Bias vectors with random numbers. 

Step 3: With the initialized values compute the output vector corresponding 
to the input vector [11 12 1з] as given below. Let the desired target 
vector be [t; 6] 





h;7funcl (wii*ii FW51* 127 W3 15 4 bn) 
017 func2 (7% tbi) 
05— func2 (wi? *h; tb;) 





Step 4: Adjustment to weights 
a) Between Hidden and Output layer. 


wy (ot+1) = wii (a) + yo (t1-01) *hi 
Wi? (n1) = wio'(n) + Yo (t2-02) *hy 








28 Chapter 1 
b) Between Input layer and Hidden layer 
Wi (0+1) = зу (n) + Yu * e* i; 


W21 (n+1) = W21 (n) T YH * e * 12 
W31 (n+1) = W31 (n) T Үн * е Т 13 








“е? is the estimated error at the hidden layer using the actual error 
computed in the output layer 


(1.е.) е- hi ж (1-h1)* [(t;-01) + (t2-02)] 
Step 5: Adjustment of Bias 
bi (n+1)= bh (n) +e* YH 


bi (п+1) i bi (n) E (ti-01)* Yo 
by (п+1) = b; (n) + (t2-02)* yo 








Step 6: Repeat step 3 to step 5 for all the pairs of input vector and the 
corresponding desired target vector one by one. 

Step 7: Let di, d2 d; ...d „бе the set of desired vectors and yı, у» уз y, be the 
corresponding output vectors computed using the latest updated 
weights and bias vectors. The sum squared error is calculated as 


SSE= (di- у1)? (do- уз)? Hds- уз)? Hda- ya)? +... (da- Yn)” 
Step 8: Repeat the steps 3 to 7 until the particular SSE is reached. 


Note if momentum is included during training updating equations are 
modified as follows. 


Wi (ntl) = wii'(n) + Yo (t1-01) * + po * Ам) (n) 


AW; (0+1) 


wi? (n1) = wi?'(n) + Awi?'(n*1) + po * Awi2(n) 


wii (ntl) = зу (п) + Awi(n+1) + pu * Aw) 
W21 (n+1) =W21 (n) nr Aw (n*l) nU Он Е Aw;i(n) 
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W31 (1+1) = w31 (п) + Aw3i(n+1) + py * Awai(n) 
bi (n+1)= bi (n) te* Үн + PH * Abi (n) 


А bi (n+1) 





b, (n 1) x b, (n) А b, (n1). po Е АБ, (n) 














b; (n 1) =, (п) А b2 (0+1): ро * АБ, (п) 
рн and po аге ће momentums used in the hidden and output layer 
respectively. 
4.3 Example 


Consider the problem for training the Back propagation Neural Network 
with Hetero associative data as mentioned below 


Input Desired Output 
000 00 
001 01 
010 01 
011 01 
100 01 
101 01 
110 01 
111 11 


ANN Specifications 


Number of layers = 3 

Number of neurons in the input layer = 3 

Number of neurons in the hidden layer =1 

Number of neurons in the output layer = 2 

Learning rate =0.01 

Transfer function used in the hidden layer = ‘logsig’ 

Transfer function used in the output layer = ‘linear’ (i.e) output is taken as 
obtained without applying non-linear function like ‘logsig’,’ tansig’. 
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Figure 1-12. Illustration of SSE decreases with Iteration 


(See Figure 1-10) 


The Figure 1-12 describes how the Sum squared error is decreasing as 
iteration increases. Note that sum squared error is gradually decreasing from 
2.7656 as iteration increases and reaches 1.0545 after 1200 Iteration. 


Trained Weight and Bias Matrix 


W1=| 1.0374 
0.9713 
0.8236 

В1= [-0.5526] 


W2= [0.8499 1.3129] 
В2- [-0.4466 -0.0163] 


OUTPUT obtained after training the neural network 
-0.1361 0.4632 
0.0356 0.7285 
0.0660 0.7756 
0.2129 1.0024 
0.0794 0.7962 
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0.2225 1.0172 
0.2426 1.0483 
0.3244 1.1747 


After rounding OUTPUT exactly matches with the Desired Output as shown 
below 
0 


OCOCOOOOCOO 
кі кі кі ка ка ка ка ОО 


4.4 M-program for Training the Artificial Neural 
Network for the Problem Proposed in the Previous 
Section 





anngv.m 

SSE=[ ]; INDEX=[ ]; 

INPUT=[0 0 0:00 1:0 1 0:0 1 1;1 00;1 0151 1 0;1 1 1]; 

DOUTPUT-([0 0:0 1:0 1 ;0 1:0 1:0 1:01 ;1 1]; 

Wi-rand(3,1); W2=rand(1,2); 

Bl=rand(1,1); B2=rand(1,2); 
for i-1:1:1200 

for j=1:1:8 

H=logsiggv(INPUT* W |+repmat(B1,[8 1])); 
OUTPUT-H*W2 +repmat(B2,[8 1]); 
ERR. HO-DOUTPUT - OUTPUT; 
ERR IH=H.*(1-H).*(sum(ERR HO')); 
Ir_ho=0.01; 
Ir_ih=0.01; 
W2(1,1)=W2(1,1)+ERR_HOG,1)*HG,1)*lr_ho; 
W2(1,2)-W2(1,2)YERR. HO(j,2)*H(j,1)*lr ho; 
WI1(1,1)2W1(1,D)*ERR IH(j,1)*INPUT(.,1)*lr ih; 
W1(2,1)=W1(2,1)+ERR IH(j,1)*INPUT(j,2)*lr ih; 
WI1(3,1)2W1(3,1)*ERR IH(j,1)*INPUT(j,3)*lr ih; 
BI-BI-ERR IH(j,1)*lr ih; 
B2(1,1)=B2(1,1)+ERR HO(j,1)*lr ho; 
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B2(1,2)=B2(1,2)+ERR HO(,2)*Ir. ho; 
end 

INDEX=[INDEX< i]; 

SSE=[SSE sum(sum(ERR HO.^2))]; 

plot(INDEX,SSE,'r'); 

xlabel(ITERATION') 

ylabel('SSE') 

pause(0.2) 

end 
H=logsigINPUT*W1+repmat(B1,[8 1])); 
OUTPUT-H*W2 +repmat(B2,[8 1]); 





logsiggv.m 


function [res]-logsiggv(x) 
res=1/(1+exp(-x)); 


5. FUZZY LOGIC SYSTEMS 


Fuzzy is the set theory in which the elements of the set are associated with 
fuzzy membership. Let us consider the set of days = (Sunday, 
Monday, Tuesday, Wednesday, Thursday, Friday, and Saturday}. This is the 
set in which elements belongs to the set with 100% membership. 

Let us consider the set of rainy days ={Sunday (0.6), Monday (0.4), 
Tuesday (0.4), Wednesday (0.6), Thursday (0.2), Friday (0.6), Saturday 
(0.3)}. This is the set in which elements are associated with fuzzy 
membership. Sunday (0.6) indicates that the Sunday belongs to the set of 
rainy days with membership value 0.6. This is different from probability 
assignment because in case of probability assignment, Sunday belongs to 
non-rainy days with membership value 0.4. (ie) (1-0.6). But in case of fuzzy 
sets Sunday belongs to non-rainy days with membership value of not 
necessarily 0.4. It can be any value between 0 to 1. 


5.1 Union and Intersection of Two Fuzzy Sets 


Let us consider the example of computing union and intersection of fuzzy 
sets as described below. 


A={a (0.3), b (0.8), c (0.1), d (0.3)) 


B- (a (0.1), b (0.6), c (0.9), e (0-.1)} 
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The elements of the set AUB are the unique elements collected from both 
the sets A and B. The membership value associated with the elements of the 
set AUB is the maximum of the corresponding membership values collected 
from the set A and B. 


Thus AUB = {a (0.3), b (0.8), c (0.9), d (0.3), e (0.1)) 


Similarly the elements of the set Ас\В are the common elements 
collected from the sets A and B. The membership value associated with the 
elements of the set AUB is the minimum of the corresponding membership 
values collected from the set A and B. 


Thus ANB= {а (0.1), b (0.6), c (0.1) 
5.2 Fuzzy Logic Systems 


Let us consider the problem of obtaining the optimum value of the 
variable ‘Z’ obtained as the output of the fuzzy logic system decided by the 
values of the variable ‘X’ and ‘Y’ taken as the input to the fuzzy logic 
system using fuzzy logic algorithm. 

The values of the variable ‘X’ vary from Xin to Xmax. Actual value of 
the variable is called crisp value. The group of the values ranges from X min to 
Xmax 1$ grouped into 3 sets XSMALL, XMEDIUM, XLARGE. The 
relationship between the crisp value and fuzzy membership value is given in 
the figure 1-13. 


X 


Fuzzy Logic 


Y System 





Figure 1-13. Fuzzy system 
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Crisp Value 


Figure 1-14. Relaionship between crisp value and Fuzzy membership value for the variable X 


The crisp values ranges from 0 to x2 belong to the set Xsmall. The crisp 
values ranges from x1 to x3 belong to the set Xmedium. The crisp value 
ranges from x2 to x4 belongs to the set Xlarge. 

The crisp value from x1 to x2 of the Xsmall set and Xmedium set with 
different membership value as mentioned in the Figure 1-14. Corresponding 
membership value decreases gradually from 1 to 0 in the Xsmall set. But the 
corresponding membership value increases gradually from 0 to 1 in the 
Xmedium set. 

Similarly the relationship between the crisp value and Fuzzy member 
ship value for the variable ‘Y’ and the variable ‘Z’ are mentioned in the 
Figure 1-15 and Figure 1-16 respectively. 
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Figure 1-15. Relaionship between crisp value and Fuzzy membership value for the variable Y 


1. Artificial Intelligence 35 





Crisp Value 


Figure 1-16. Relaionship between crisp value and Fuzzy membership value for the variable Z 


The identical relationship between the crisp value and fuzzy value of all 
the variables X, Y and Z are chosen for the sake of simplicity (see figure). 
Steps involved in obtaining the crisp value of the output variable for the 
particular crisp values for the variable X and Y are summarized below. 


5.2.1 Algorithm 


Step 1: Crisp to fuzzy conversion for the input variable ‘X’ and ‘Y’ 


Let us consider the crisp value of the input variable X be Xinput the crisp 
value of the input variable Y be Yinput. Let Xinput belongs to the set 
Xsmall with fuzzy membership 0.7 (say) belongs to set Xmedium 0.3 (say) 
as shown in the Figure 1-17 and Figure 1-18 given below. 
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Figure 1-17. Xinput as the crisp input for the variable X 
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y4 
Crisp Value 


Figure 1-18. Yinput as the crisp input for the variable Y 


Similarly the Yinput belongs to the set Ymedium with fuzzy membership 0.7 
and belongs to the set Ylarge with 0.2 as shown in the figure given below 


Step 2: Fuzzy input to fuzzy output mapping using fuzzy base rule 


Fuzzy base rule is the set of rules relating the fuzzy input and fuzzy output 
of the fuzzy based system. Let us consider the Fuzzy Base Rule as shown in 
the table given below. 


Table 1-1. Fuzzy Base Rule 








Ysmall Ymedium Ylarge 

small Zmedum Zmedium Zsmall 

Xmedium Zmedium Zlarge Zlarge 
Xlarge Zlarge Zlarge Zmedium 





Thus Fuzzy input set Xsmall = {Xinput (0.7) } 
Xmedium = {Xinput (0.3) } 
Xlarge = {Xlarge( 0 ) 


Ysmall = {Yinput (0) } 

Ymedium = {Yinput (0.7) } 

Ylarge = {Xlarge(0.2) 
Fuzzy output set is obtained in two stages. 


Intersection of Fuzzy sets 

Xsmall ^ Ysmall = Zmedium ={Zoutput(0)} 
Xsmall ^ Ymedium = Zmedium ={Zoutput(0.7)} 
Xsmall ^» Ylarge = Zsmall ={Zoutput(0.2)} 
Xmedium r^ Ysmall = Zmedium ={Zoutput(0)} 
Xmedium r^ Ymedium = Zlarge={Zoutput(0.3)} 


e о 6 — 
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e Xmedium n Ylarge = Zlarge ={Zoutput(0.2)} 

e Xlargem Ysmall = Zlarge ={Zoutput(0)} 

e Xlargern Ymedium = Zmedium ={Zoutput(0)} 

e Xlargem Ylarge = Zsmall ={Zoutput(0)} 

2. Union of Fuzzy sets 


Zsmall = U(Zoutput(0.2), Zoutput(0))={ Zoutput(0.2) } 
Zmedium = U(Zoutput(0), Zoutput(0.7), Zoutput(0),Zoutput(0)) 
= {Zoutput (0.7)} 
e Zlarge = U(Zoutput(0.3), Zoutput(0.2), Zoutput(0))={Zoutput(0.3)} 


Thus Fuzzy output based on the framed fuzzy base rule is given as 


Zsmall = {Zoutput(0.2) 
Zmedium = {Zoutput(0.7} 
Zlarge = {Zoutput(0.3)} 


Step 3: Fuzzy to Crisp conversion for the output variable Z using 
Centroid method 
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Figure 1-19. Defuzzification 
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Truncate the portions of the Zsmall, Zmedium and Zlarge sections of the 
output relationship above the corresponding thresholds obtained from output 
fuzzy membership as shown in the figure. Remaining portions are combined 
with overlapping of the individual sections. Combined sections can also be 
obtained by adding the individual sections instantaneously as illustrated in 
the example given in the section 5.4 Obtain the crisp value Zout which 
divides the combined sections into two halves as shown in the figure. Thus 
crisp value Zout is obtained for the corresponding input crisp value Xinput 
and Yinput using Fuzzy logic system. 


5.3 Why Fuzzy Logic System? 


The Automatic system uses set of rules framed with the input data to 
decide the output data. For example the automatic system is designed to 
maintain the temperature of the water heater by adjusting the heat knob 
based on the water level at that particular instant and temperature at that 
particular instant. In this example water level measured using the sensor at 
that particular instant and temperature measured using the sensor at that 
particular instant are the input to the control system. Control system uses the 
set of rules framed as given below for deciding the position of the heat knob 
which is the output of the control system 


e If T1<=Temperature <=T2 and W1<=Water level<=W2, Adjust the heat 
knob position to K1. (Say). 

е If TO<=Temperature <=T1 and W0<=Water level<=W1, Adjust the heat 
knob position to K2. (Say) etc. 


Suppose if the temperature is just less than Т1, this rule is not selected. 
The small variation ‘AT’ makes the automatic control system not to select 
this particular rule. To overcome this problem instead of using crisp value 
for framing the rule, fuzzy values can be used to frame the rules (i.e) without 
using the actual values obtained through input sensors This is called fuzzy 
base rule as given below 


e Ifthe temperature of the water is small and water level is high, place the 
heat knob in the high position. (say) 

e If the temperature of the water is medium and water level is less, place 
the heat knob in the less position. (say) etc. 
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In this case the problem only fuzzy sets are used for decision making which 
makes the system to take decision like the human decision. 

To design the fuzzy based system, the knowledge about the system 
behavior is reguired to find the relationship between crisp value and the 
fuzzy value (see above) and to frame the fuzzy base rule. This is the 
drawback of the fuzzy based system when compared to ANN based system 
which doesn’t require any knowledge about the system behavior. 


5.4 Example 


Let us consider the problem of obtaining the optimum value of the 
variable ‘Z’ obtained as the output of the fuzzy logic system decided by the 
values of the variable ‘X’ and ‘Y’ taken as the input to the fuzzy logic 
system using fuzzy logic algorithm. The crisp -fuzzy relationship is as 
mentioned in the section 5.2. The fuzzy base rule is as given in the table 1-1. 

The values for the variables used in the section 5.2 are assumed as 
mentioned in the table 1-2. 


Table 1-2. Assumed values for the variables used in the section 5.2 





Xl х2 X3 х4 Ү1 Y2 Y3 Y4 Zl Z2 Z3 Z4 
3 5 7 10 20 45 65 100 35 70 90 100 








For the arbitrary input crisp values xinput=6 and yinput=50, zout is obtained 
as 78.8 using fuzzy logic system described in the section 5.2 


Xinput 

















0 1 2 3 4 5 6 7 8 9 10 
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Figure 1-20. Relationship between crisp value and fuzzy membership for the Input variable X 
in the example 
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Figurel-21. Relationship between crisp value and fuzzy membership for the Input variable Y 
in the example 
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Figure 1-22. Relationship between crisp value and fuzzy membership for the Input variable Z 
in the example 
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Figure 1-23. Defuzzification in the example 


1. Artificial Intelligence 41 


5.5 M-program for the Realization of Fuzzy Logic 
System for the Specifications given in Section 5.4 





fuzzygv . m 


posx=[0 3 5 7 10]; 

maxx=max(posx); 

posx=posx/max(posx); 

posy=[0 20 45 65 100]; 

maxy=max(posy); 

posy=posy/max(posy); 

posz=[0 35 70 90 100]; 

maxz=max(posz); 

posz=posz/max(posz); 

i=posx(1):1/999:posx(2); 

J=(posx(2)):(1/999):posx(3); 

k=(posx(3)):(1/999):posx(4); 

l=(posx(4)):(1/999):posx(5); 

xsmall=[1/posx(2)*i (-1/(posx(3)-posx(2))*(j-posx(3))) zeros(1,length([k 1]))]; 

xmedium =([zeros(1,length([i])) (1/(posx(3)-posx(2))*(j-posx(2))) (-1/(роѕх(4)- 
posx(3))*(k-posx(4))) zeros(1,length([1]))]; 

xlarge=[zeros(1,length([i j])) (1/(posx(4)-posx(3))*(k-posx(3))) (-1/(роѕх(5)-роѕх(4))*(1- 
posx(5))) ]; 


figure 

plot([i j k 1]*maxx,xsmall,':'); 
hold 

plot([i j k 1] *maxx,xmedium,-"); 
plot([i j k 1]*maxx,xlarge,'--'); 
xlabel('crisp value") 
ylabel('fuzzy value") 

title" XINPUT") 


i=posy(1):1/999:posy(2); 

j=(posy(2)):(1/999):posy(3); 

k=(posy(3)):(1/999):posy(4); 

l=(posy(4)):(1/999):posy(5); 

ysmall=[1/posy(2)*1 (-1/(posy(3)-posy(2))*(j-posy(3))) zeros(1 ,length([k 1]))]; 
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ymedium =[zeros(1,length([i])) (1/(роѕу(3)-роѕу(2))*()-роѕу(2))) (-1/(роѕу(4)- 
роѕу(3))*(К-роѕу(4))) zeros(1.length([1]))]; 

ylarge=[zeros(1 ,length([i j])) (1/(posy(4)-posy(3))*(k-posy(3))) (-1/(роѕу(5)-роѕу(4))*(1- 
posy(5))) 1; 


figure 

plot([i j k 1]*maxy,ysmall,':'); 
hold 

plot([i j К 1]*maxy,ymedium,'-); 
plot([i j k 1]*maxy,ylarge,'--'); 
xlabel('crisp value") 
ylabel('fuzzy value") 

title" YINPUT") 


i=posz(1):1/999:posz(2); 

J=(posz(2)):(1/999):posz(3); 

k=(posz(3)):(1/999):posz(4); 

l=(posz(4)):(1/999):posz(5); 

zsmall=[ 1/posz(2)*i (-1/(posz(3)-posz(2))*(j-posz(3))) zeros(1,length([k 1]))]; 

zmedium =[zeros(1,length([1])) (1/(posz(3)-posz(2))* G-posz(2))) (-1/(posz(4)- 
posz(3))*(k-posz(4))) zeros(1,length([1]))]; 

zlarge=[zeros(1,length([i j])) (1/(posz(4)-posz(3))*(k-posz(3))) (-1/(роѕ2(5)-роѕ2(4))*(1- 
posz(5))) ]; 


figure 

plot([i j k 1] *maxz,zsmall,':"); 
hold 

plot([i j k 1]]žmaxz,zmedium,"-"); 
plot([i j k 1]*maxz,zlarge,'--'); 
xlabel('crisp value') 
ylabel('fuzzy value") 

title" ZOUTPUT) 


xinput=[6]; 
xinput=round(((xinput/maxx)*1000)) 
yinput=[50]; 
yinput=round(((yinput/maxy)* 1000)) 
Vofuzzy values 


xfuzzy=[xsmall(xinput) xmedium(xinput) xlarge(xinput)]; 
yfuzzy=[ysmall(yinput) ymedium(yinput) ylarge(yinput)]; 
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%Output fuzzy values obtained using framed fuzzy base rule 

zfuzzys-max([(min([xfuzzy(1) yfuzzy(3)]))]); 

zfuzzym=max([(min([xfuzzy(1) yfuzzy(2)])) (min([xfuzzy(2) yfuzzy(1)])) 
(min([xfuzzy(3) yfuzzy(3)]))]; 

zfuzzyl=max([(min([xfuzzy(2) yfuzzy(2)])) (min([xfuzzy(2) yfuzzy(3)])) (min([xfuzzy(3) 
yfuzzy(1)])) (min([xfuzzy(3) yfuzzy(2)]))]); 


%Defuzzification 

zs=reshape((guantiz(zsmall,zfuzzys)*2-2)*(- 
1/2),1,size(zsmall,2)).*zsmall+quantiz(zsmall,zfuzzys)'.*zfuzzys; 

zm-reshape((quantiz(zmedium,zfuzzym)*2-2)*(- 
1/2),1,size(zmedium,2)).*zmedium--quantiz(zmedium,zfuzzym)'.*zfuzzym; 

zl=reshape((quantiz(zlarge,zfuzzyl)*2-2)*(- 
1/2),1,size(zlarge,2)).*zlarge+guantiz(zlarge,zfuzzyl)'. *zfuzzyl; 


figure 

plot([i j k 1]*maxz,zs,':') 
hold 

plot([i j k 1]*maxz,zm,'-') 
plot(zl,'--') 


final-sum([zs;zm;zl]); 


figure 
plot([i j k 1]*maxz ,final) 
title DEFUZZIFICATION") 


S=cumsum(final) 
[p.q]=find(S>=(sum(final)/2)) 
zout=(q(1)/1000)*maxz; 

hold 
plot(ones(1,100)*zout,[0:1/99:1],"*") 
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6. ANT COLONY OPTIMIZATION 


Ant colony optimization technigues exploit the natural behavior of ants in 
finding the shortest path to reach destination from the source. 

The figure 1-24 is the top view of the ants moving from source to 
destination in search of food. There is the obstacle in the middle. There are 
two paths available so that the ants can choose to move from source to 
destination. As expected after some time duration ants prefer to choose 
path1, which is the shortest distance between source and destination. This is 
due to the fact that ants excrete the pheromones when they are moving. As 
more ants have crossed the shortest path within the stipulated time when 
compared to the other path, more pheromones are excreted in the shortest 
path. This helps the following ants to choose the shortest path to cross the 
obstacles. 

The mathematical eguivalent of this behavior is used in optimization 
technigue as the Ant colony optimization. 


6.1 Algorithm 


Consider the problem of finding the optimum order in which the numbers 
from | to 8 are arranged so that the cost of that order is maximized. Cost of 
the particular order is computed using the two matrices A and B as described 
below. 





Figure 1-24. Illustration of Ant Colony 
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Figure 1-25. Matrix A 
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Figure 1-26. Matrix B 


The cost ofthe order [32157246 8] is computed as the sum of the 
product of values obtained from the positions (1,3) or (3,1), (2,2), (3,1), (4,5) 
or (5,4) , (5,7) or (7,5) ,(6,2), (7,4) and (8,6) of the matrices A and B . 

Cost ([32 1572 46 8] ) = a6*b6 + a3*b3 + a6*a6 + a14*b14 + a26*b26 + 
a17*b17 +a25*b25 + a36*b36 
(see Matrix A and Matrix B ) 


The problem of finding the optimal order is achieved using Ant colony 
optimization technique as described below. 


Step 1: Generate the 5 sets of orders randomly which are treated as the 
initial values selected by the 5 Ants respectively. 


The following are the orders selected by the 5 ants along with the 
corresponding cost measured as described in the previous section. 
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Table 1-3. Initialization of 5 values to the 5 Ants (say) 

ANT ORDERS COST 

ANT 1 71362548 Cl 

ANT 2 24863751 C2 

ANT 3 723 516 C3 

ANT 4 15487326 C4 

ANT 5 48536271 C5 





Step 2: Pheromone is the value assigned for selecting the i" number for the 
j^ position of the sequence selected by all the ants.(i.e) it is the 
matrix in which the value at the index (i ,j) is the value assigned for 
selecting the ith number for the jth position of the order selected by 
all the ants. Initially the values of the matrix are completely filled up 
with ones. The values of the Pheromone matrix are updated in every 


iteration as 


Pheromone matrix (n+1) = Pheromone matrix (n) + Updating Pheromone matrix 


The value of the updating pheromone matrix must be high if the minimum 
cost is assigned. Hence the matrix is filled up with the inverse values of the 
corresponding cost as displayed below. 


Table 1-4. Updating Pheromone Matrix 








1 2 3 4 5 6 7 8 
1 1/C4 1/С1 0 0 1/C3 0 0 1/C2+1/C5 
2 1/C2 1/C3 0 0 1/С1 1/С5 1/С4 0 
3 0 0 1/C1+1/C3 1/С5 1/C2 1/C4 0 0 
4 1/C5 1/C2 1/С4 0 0 0 1/C1+1/C3 0 
5 0 1/С4 1/С5 1/C3 0 1/С1 1/C2 0 
6 0 0 0 1/C1+1/C2 1/С5 1/C3 0 1/C4 
7 1/C1+1/C3 0 0 0 1/С4 1/C2 1/C5 0 
8 0 1/С5 1/C2 1/C4 0 0 0 1/C1+1/C3 
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Step 3: Next set of orders selected by the ants is as described below. 


e Generate the position sequence randomly 
7 8 2 3 6 4 1 5(say) 


Let the order decided by the antl be represented as [ul u2 u3 u4 u5 иб 
17118 | 


e The value for the u7 in the order as mentioned above is selected first as 
the first number of the randomly generated position sequence is 7. 

e To select the value for the 17, 7" column of the pheromone matrix is 
used .(Note that 7^ column belongs to 7" position of the order generated) 

e 7° column of the pheromone matrix is rewritten below as Column7 
vector arranged in row wise for better clarity. 


Column7 -[0 1/C4 0  I/CI-1/C3 1/C2 0 1/С5 0] 


е Normalized value of the 7" column of the pheromone matrix is given 
below as Normalized vector7 arranged in row wise. This is also called as 
the probability vector used for selecting the value for the u7 in the order. 


S = [0 + 1/C4 + (1/C1+1/C3) + 1/C2 + 1/C5] 
Normalized vector7 = (Column 7) / S = [pl p2 p3 p4 p5 p6 p7 p8] (say) 





Select the number corresponding to the position of the lowest value of the 
normalized vector. Suppose if p5 is the lowest among the values in the 
vector, the number 5 is assigned to the variable u7. 


e The next value in the position sequence is 8.Hence the value for the 
variable u8 15 selected. But the value cannot be 5 as it was already 
assigned to the variable u7. 8^ column of the pheromone matrix except 
the 5" row is used for selecting the value to be assigned to the variable u8 
as described below. 

e 8^ column of the pheromone matrix is rewritten below as Column 8 
vector with 5" row filled up X indicating 5" row is not consider for 
selection. The vector is rearranged in row wise as displayed below 


Column 8= [(1/C2 +1/С5) 0 00 X 1/C4 0 (I/CI-I/C3)] 


е Normalized value of the 8" column of the pheromone matrix is 





S-[(1/C2 +1/C5)+ 0 +0 +0 +1/C4+0+(1/C1+1/C3)] 


Normalized vector8 = (Column 8) / S = [q1 q2 43 q4 q6 97 98] 
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Select the number corresponding to the position of the lowest value of the 
Normalized vector8. Suppose if 43 is the lowest among the values in the 
vector, the number 3 is assigned to the variable u8. 

This process is repeated to obtain the orders selected by the Ant 2, Ant3, 
An4 and Ant 5.This is called one iteration. 


Step 4: 

e The updating pheromone matrix for the second iteration is computed as 
described in the step 2 using the set of orders selected by the five ants in 
the first iteration. This is called as Updating pheromone matrix (2). 

е Pheromone matrix used in the 2" iteration is computed as 


Pheromone matrix (2) = Pheromone matrix (1) +Updating Pheromone matrix (2) 


Step 5: Next set of orders selected by the five ants are computed as 
described in the step 3. 


Step 6: The step 3 to step 5 is repeated for the particular number of 
iterations. Thus best set of orders selected by the five ants are 
obtained using Ant colony optimization. 


Note that best among the best set of orders selected by the five ants 
in the last iteration is selected as the global best order. 


6.2 Example 


Consider the problem of finding the optimum order in which the numbers from 
1 to 8 are arranged so that the cost of the order is maximized. Cost of the 
particular order is computed as described in the section 6.1. The matrices A and 
B are given below. 


A= 
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Initial orders selected by the four ants are displayed below. 
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Number of ants used to find the optimum order is 4. The orders selected 


by the four ants after 50 Iterations are displayed below. 


ANTI: 
ANT2: 
ANT3: 
ANTA: 


1 
1 
3 
4 
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92-4 92 99 


RR 


— к. Un л 
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00 OO OO OO 


Note that ANTI and ANT2 found the best order as expected using Ant 


colony technique. 


Initially the cost values are having more changes. After 40" iteration, 
cost value gradually increases and reaches maximum which is the optimum 
cost corresponding the optimum order selected using Ant colony 


technique. 
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Figure 1-27. Iteration (vs.) Best Cost selected by the 4 ants 


6.3 M-program for Finding the Optimal Order using Ant 
Colony Technique for the Specifications given in the 
Section 6.2 





antcolonygv.m 


% Consider the problem of finding the optimum order in which the numbers from 1 to 8 
are arranged so that 

Хоће cost of that order is decreased .Cost of the particular order is computed using the 
two matrices A and B as described below. 

%The cost of the order [32157246 8] is computed as the sum of the product of 
values 

%obtained from the positions (1,3) or (3,1),(2,2),(3,1) ,(4,5) or (5,4) , (5,7) or (7,5) ,(6,2), 
(7,4) and (8,6) of the matrices A and B . 

%Cost ([3 21572468] ) =a6*b6 + a3*b3 + a6*a6 + a14*b14 + a26*b26 + al7*b17 
+a25*b25 + a36*b36 see Matrix A and Matrix B 








%Form the matrices given below 
matA-diag([ones(1,9)*9]); 
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matB=diag([ones(1,9)*9]); 
%Initializing the orders selected by 4 Ants 


[p1, ANT 1 J=sort(rand(1,8)); 
[p2, ANT2]=sort(rand(1,8)); 
[p3, ANT3 J=sort(rand(1,8)); 
[p4, ANT4]=sort(rand(1,8)); 
Y%Intializing the Pheromone matrix 
% columns are the positions 
Vorows are the value of the position 
ANTS=[ANT1;ANT2;ANT3;ANT4]; 
P=ones(8,8); 
col=[]; 
figure 
hold 
for 1=1:1:50 
%Updating pheromones 
cl=computecost(ANTS(1,:),matA,matB); 
c2=computecost(ANTS(2,:),matA,matB); 
c3=computecost(ANTS(3,:),matA,matB); 
c4=computecost(ANTS(4,:),matA,matB); 
col=[col max([c1 c2 c3 c4])]; 
plot(col) 
pause(0.01) 
c=[cl c2 c3 с4]; 
for 1=1:1:8 
for j=1:1:8 
[х,у] = find(ANTS(,1)==]); 
for k=1:1:length(x) 
Р(,)=Р(,1)-+(1/е(х(Ю))); 
епа 
епа 
епа 


АМТЅ=[]; 
for 1=1:1:4 
P cur=P; 
tempant=zeros(1,8); 
notation=[1 234567 8]; 
[p.position]=sort(rand(1,8)); 
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for 1=1:1:8 

temp-P cur(:,position(1)); 

temp=temp/sum(temp); 

[value,index J=min(temp); 

tempant(position(i))=notation(index); 

P_cur(index,:)=[]; 

notation(index)=[]; 

end 

ANTS=[ANTS;tempant]; 

end 

cl=computecost(ANTS(1,:),matA,matB); 
c2=computecost(ANTS(2,:),matA,matB); 
c3=computecost(ANTS(3,:),matA,matB); 
c4=computecost(ANTS(4,:),matA,matB); 
end 





computecost.m 
function [s]=computecost(A,P,Q) 
s=0; 
for 1=1:1:8 
s=s+P(max(i,A(i)),min(i,A(i))).*Q(max(i,A(i)),min(i,A(i))); 
end 
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PROBABILITY AND RANDOM PROCESS 
Algorithm Collections 


1. INDEPENDENT COMPONENT ANALYSIS 


Consider “т” mixed signals yi(t), ya(t)...Ym(t) obtained from the linear 
combination of ‘n’ independent signals x,(t), x(t) ... x,(t), where n>=m is 
given below. 


yi(t) = 811 xi(t) т 412 х›(0) T 413 xa(t) pp .атХха(0 
y2(t) = a21 X(t) т 422 х›(0) т 423 x(t) т ...аљһхХь (t) 
y3(t) — 831 xi(t) T 432 х›(0) т 433 xa(t) Run .a3nXn(t) 
Valt) = ац Xi(t) + ад X2(t) + адз X3(t) + .. .a4nX (0) 











Ym(t) = amı Xi(t) + am KAL) + Ama X3(t) + ...AmnXrt) 


Independent signals are obtained from the mixed signals using 
Independent Component Analysis (ICA) algorithm as described below 


1.1 ICA for Two Mixed Signals 


The samples of the signals xi(t) and x2(t) are represented in the vector 
form as ху=[ху Xi? X13 ...Xi4] and x?7[xoi X2» X23 ... X2n] respectively. The 
signals yı(t)=a1ı*xı(t)+aı12*x2(t) and y2(t)=a21*xi(t)+a22*x2(t) are represented 
in the vector form as У1=[у1 yi2 Уіз Yin] and yo=[V21 yo Yz ... Уа] 
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respectively. Thus the ICA problem сап be represented in the form of matrix 
as described below. 


yu Уз! X11 X21 

Уі yo X12 X22 

yis yos X13 X23 

Yi4 y24 X14 X24 а аш 
yis yas Х15 X25 

У16 Y26 = X16 X26 а: ао 
Yin Уһ Xin X2n 

[Y] = [X] [A] 


Given the matrix [Y], obtaining the matrices [X] and [A] such that the 
columns of the matrix [Х] are independent to each other is the ICA problem. 

Consider three independent signals and the corresponding mixed signals 
along with the kurtosis values are displayed below. 


15 1 
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Kurtosis is the statistical parameter used to measure the Gaussian nature 
of the signal. Kurtosis is inversely proportional to the Gaussianity of the 
signal. Note that the kurtosis values are maximum for independent signals 
compared to the mixed signals. (i.e.) Independent signals are more non- 
Gaussian compared with mixed signals. Mathematically kurtosis is 
computed using the formula as displayed below, where E[X] is the 
expectation of the vector X 


ICA problem is to obtain the matrices [X] and [A] such that column 
vectors of the matrix [Х] are independent to each other. (i.e.) Kurtosis 
values computed for the column vectors are maximum. 


[Y] = IX] [AT 
(ie) [Y] = [A] [X] 
ІХІ-ТАТ [Y] 
[X]-EB]LY] 
XI and X2 are the two column vectors of the matrix [X]' 
(i.e) [X] = [X1 X2] 
Kurtosis measured for the Column vector [X1] is given as 


E [X1*}-3{E [X]? ? 


Similarly kurtosis measured for the column vector [X2] is given as 


Е [X2^-3(E [X27]? 
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Estimating the best values for the matrix [В] such that kurtosis measured 
for the column vectors of the matrix [X]! = [[B][Y]]' as defined above is 
maximum 

We require another constraint about the matrix [B] to obtain the values of 
the matrix [B] which is described below 

The covariance matrix (COV(Y)) computed for the group of vectors 
collected row-by-row of the matrix Y' is computed using the formula as 
displayed below 





Let M = | (yu + y12 tyist...Yin) /n 
(Yor + yz +V23t+...Y2n)/n 


COV(Y) = 
Уп yu Ул | +y ||Уо Yo} +..| Уш Yin Y2n|| m - ММ" 
y21 yo Yon 


(іе) COV(Y) = E[Y Y7] - MM" 


The matrix computed is of size 2x2. The value at the position (1,1) 
conveys the information about how the first elements of the collected vectors 
are correlated with itself (variance). The value at the position (1,2) conveys 
the information about how the first elements are correlated with second 
elements of the collected vectors. 

For example the covariance matrix computed for 2D vectors collected 
from the particular independent signals and the corresponding mixed signals 
are given below. 

0.2641 x 10° — -0.1086 x 10° 

COV independent = 

-0.1086 x 10° 0.5908 х 107% 


COV mied = 0.1086 0.0613 


0.0613 0.0512 
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To compare the covariance matrix, correlation coefficient is used. It is the 
normalized covariance matrix whose (m, n) element is computed as 
COV( m, n) / (sqrt(COV(m,m)*COV(n,n)) 

The correlation co-efficient for the above covariance matrix (COV) is given 
below 


CORR independent = 1 -0.2750 
-0.2750 1 


CORR mixed = |l 0.822 


0.822 1 


Note that cross correlation value is less for independent signals compared to 
the mixed signals. This fact is used as the source for second constraint. (1.е.) 
the correlation coefficients matrix of the independent signals is almost 
identity matrix. Also if the variances of the signals are unity the covariance 
matrix and the correlation co-effients are identical. 

The mean and variance of the mixed signals are made 0 and 1 respectively 
so that the mean and variance of the row vectors of the matrix [Y] are 0 and 
1 respectively. [Refer Mean and variance Normalization Section 4 of 
Chapter 2] 

The matrix [Y] can be transformed into another matrix [Z] such that the 
covariance matrix computed for the 2D vectors collected from the matrix [Z] 
as described above is almost unit vector (diagonal) using KLT (Refer 
Hotelling transformation). (i.e.) E[[Z][Z] ]H[1] 

Let us define the transformation matrix [T], which transforms the matrix 
[Y] into [Z] as described below. 


[Z]=[T][Y] implies 
[Y]H[T] ' [Z] implies 
[Y ]-[U] [Z] 


Also, [X]-[B][ Y] 
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[X]-[B][U][Z] 


[B][U][Z]=[X] 
[BIUIIBIIB]'[ZJ=[X] 


[A][Z]={[B][U][B]}" [X] 


Note that A=[B]' 


Column vectors of the transformation matrix [U] are orthonormal to 
each other. (ie) [U] ' = [U]' 


If E [XX'] = [I] (Identity Matrix) 
((1.e.) Covariance matrix of the Independent signals with mean zero and 
variance one is Identity matrix ) 


Computing Covariance Matrix of the RHS of the equation 
[A] [Z]={[B][U] [B]}" IX] 


ЕІСІВІІЛІВІУ [XD СІВІЛІВІ; XD] 


= E [(«[B]TU]JIBI D [X] EX] СІВІСІВІ >" ] 


= ({[B][U][B]}") E [XX'] (£([B]TU][BIY )" 


= (([BJ[UI[B]3”) ({[в]ГО[в] >" 


= (BT [UT [BT CIBT [UT [BT )' 


-[BT [UT [BT [E [U ] "B" 
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If the columns of the В matrix is orthonormal to each other 
(i.e) [B] ' = [B] 


[B] [UT [BI [B] [U ] [в] 


-[BT [UT [U ] "By" 


= [B] [BJ (Because [U]'=[U]') 


-[B] [B] 


— [I] (Identity Matrix ) 


Computing Covariance Matrix of the LHS of the equation 
[A][Z]={[B][U][B]}" [X] 
E ([A][Z)) (А21) ) 
=AE[ZZ'JA" 
= [А]ЦГА]' 
-ІВГІВ T 
= Ш 
(Because [B] '=[B]' is assumed ) 
Thus Covariance matrix ofthe LHS and RHS are Identity matrix if 
е E[XX'] is the Identity matrix (Assumed) 
е E[ZZ’] is the Identity matrix (Property) 


e [B]'- [B]! (Assumed) 
e [U]'=[U]' (Property) 
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In other words, If the Columns of the matrix [В] are orthonormal to each 


other (ie) [В]! = [B]', then E[XX'J=[I] (i.e) Covariance Matrix of the 
independent signals is the Identity Matrix.(Required). 


[A][Z]={[B][U][B]}" [X] 


When [A] is estimated such that covariance matrix of the LHS and RHS are 
Identity Matrix, [A][Z]=[X] (ie) ([B][U][B]] ! becomes the Identity matrix. 


Thus ICA Problem is the optimization problem as mentioned below. 

Given the matrix [Y], estimating the best values for the matrix [B] such 
that kurtosis measured for the row vectors of the matrix [X] = 
[A][Z]=[BĦ][Z] is maximum. 


Subject to the constraint B'B=[I] (i.e) column vectors of the matrix В is 
orthonormal to each other. 


Let [B]=|b1 bi; 
ba b» 


Z= | 71 212 713... Zin 


2721 722 723... Zn 


T ad 
[B] [Z] = bi; bz 211 712 713... Zin 
bi? b» 721 722 723... Z2n 


=|buziutb2iZ21 бз1712+021722 ... biiZintb21Zon 











biZiitb22Z21 bi2Z12+b22222 ... 61221602222 


Kurtosis of the first row = E[(biiziitboizoi )] — 3 1  (Әилыы2а Y]Y? 


Kurtosis of the second row= E[(bi2z1+b2222; )“] - 3 {E[(bi2z1 о? ) 1}? 
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Constraint [B]'[B]=[I 1 


(i.e ) bi; by; +b; by) =1 
biz 6+6» =1 





Thus Iterative approach to estimate the best values for В matrix using 
Lagrange’s method is described below. 


Partial differentiating the function 

E[(biiziboizoi )] — 3 (E[(buiziitboizoi УР +A (bii birtbai b21 —1) 
with respect to the unknown constants Бі, b 21 and equate to zero. 

With respect to b;; 

E[4(b i iziitboizoi y zii ]-3 *214E[(biiziitboizoi YhE[biiziitboizo; )zi] 

TA*2*b;; =0 
E[(biiziitboizoi y] is the variance of the independent signal. This can be 
assumed to be 1. 

Also E (zii 211) is the variance of the first mixed signal and hence it is 1. E 
(zii 221) is the covariance between the two mixed signals and the value is 
Zero. 

Therefore the above equation becomes 


4*E[(biizitbo»izoi y zi ]-6*2 *b,,+A*2*b,, =0 


Let the Lagrange multiplier A be (-2). [For More Details Refer 
Estimation Techniques] 


Thus Iteration equation becomes 

bii(tt1) = E[(bii (t)zitboi(t) Z2i) zii]-3*bi (0) 
bii(t) is the value for bı; in t? iteration 

bii(t+1) is the value for b;; in (t-1)^ iteration 


Similarly iteration equation for other elements of the matrix B are 
computed and are displayed below 
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bii(t+1) => E[(bii (9711-5210 721 y Zi]-3*bii (t) 


bo (t+1) = E[(bii (Ozirtbor(t) Zz) 241-351 (t) 


b,2(t+1) Е Eb (9711-5200) 721 y 211 1-3 *bp (t) 








601) = Eb (971-8220 221 y 221 ]-3*b» (t) 


Recall that the columns of the matrix B are made orthonormal to each 
other. This can be explicitly performed in every iteration, after updating the 
values using the equation given above. 


1.1.1 


Step 1: 


Step 2: 


Step 3: 


Step 4: 


Step 5: 


Step 6: 


ICA algorithm 
Given mixed signals y1(t) and y2(t) are converted into z1(t) and 
Z2(t) using Hotellilng transformation such that the covariance matrix 
computed using the converted signals z1(t) and z2(t) is the identity 
matrix.[Use Hotelling Transformation] 


Initialize the values for the matrix B such that B'B=[I] 


Update the elements of the matrix B using the Iteration formula 
displayed above. 


Columns of the matrix B is made orthonormal to each other as 
mentioned below. 


В = B * real(inv(B' *B)*(1/2)); 
Repeat step 3 and step 4 for ‘N’ Iteration. 


Independent signals are obtained by multiplying В! with the Z 
matrix 
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Example: Independent signals x1 (t) and x2 (t) are mixed to get the 
signals y1 (t)=0.3*x1 (t)+0.7*x2 (t), y2 (t)=0.7*x1 (t)+0.3*x2 (t) as shown 
in the figure given below 


х1 (t) normalilzed to zero mean and unit variance 
2 T T T T T T T T T 











П 1 1 1 1 1 1 1 1 
ü 1000 2000 35000 4000 5000 6000 ғооо 8000 5000 10000 


x2 (0 normalized ta zero mean and unit variance 
1 T T T T T T T 














1 1 1 1 1 1 1 1 1 
ü 1000 2000 3000 4000 5000 6000 7000 8000 5000 10000 


y 10-0. 7*1 0H. 3*x2(t) and normalized to zero mean and unit variance 








1 1 1 1 1 1 1 1 1 
ü 1000 2000 5000 4000 5000 6000 “000 8000 s000 10000 


y2 (-0.37x1()-À1.7*x2(t) and normalized to zero mean and unit variance 








П П П П П П П П П 
ü 1000 2000 3000 4000 5000 6000 уйі 8000 s000 10000 
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Using ICA independent signals are obtained in two stages with 
decorrelated signals zl(t) and z2(t) in the first stage and the Independent 
estimated signals x1(t) and x2(t) in the second stage as displayed below. 
Number of Iterations (N) used 15 10000. 


zit) normalized to zero mean and unit variance 














1 1 1 1 1 1 1 
ü 1000 2000 5000 4000 5000 6000 F000 8000 3000 10000 


z2(t) normalized ta zero mean and unit variance 
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Estimated x 1) using ICA 
2 T T T T T T 
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Estimated x20) using ICA 
T T 








1 1 1 1 1 1 1 1 1 
ü 1000 2000 5000 4000 5000 6000 7000 000 s000 10000 
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1.2 M-file for Independent Component Analysis 





icagv.m 


%ICA WITH THREE COMPONENTS 


close all 

clear all 

1=0:1/10:1000; 
a=sin(0.1*pi*i); 
a=a(1:1:10000); 
a=meanvarnorm(a,0, 1) 
b=randn(1,10000); 
b=meanvarnorm(b,0,1); 
1=0:1/10:1000; 
c=exp(0.001 *pi*i); 
c=meanvamorm(c,0, 1); 
c=c(1:1:10000); 
s1=0.9*a+0.6*b+0.6*c; 
sl=meanvarnorm(s1,0,1); 
s2=0.6*a+0.9*b+0.6*c; 
s2=meanvarnorm(s2,0, 1); 
s3=0.6*a+0.6*b+0.9*c; 


s3=meanvarnorm(s3,0, 1); 


























figure 
subplot(3,1,1) 
plot(a); 
subplot(3,1,2) 
plot(b); 
subplot(3,1,3) 
plot(c); 

figure 
subplot(3,1,1) 
plot(s1); 
subplot(3,1,2) 
plot(s2); 
subplot(3,1,3) 
plot(s3); 

%ICA START 
datal=[s1;s2;s3]; 
MEANMAT=repmat(mean(data1')',1,length(data1)); 
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cl=cov(datal'); 
[T1,E1J=eig(c1); 
data2=T 1'*(datal-MEANMAT); 
%data2=T 1 *(data1); 
c2=cov(data2'); 
[T2,E2]=eig(c2); 
data-E2^(-1/2)*data2; 
T-E2^(-1/2)*TI; 
c=cov(data'); 
figure 
subplot(3,1,1) 
plot(data(1,:)); 
subplot(3,1,2) 
plot(data(2,:)); 
subplot(3,1,3) 
plot(data(3,:)); 
%ICA STARTS 
w1 l=rand(1); 
w12=rand(1); 
w13=rand(1); 
w21-rand(1); 
w22-rand(1); 
w23-rand(1); 
w31-rand(1); 
w32-rand(1); 
w33-rand(1); 
Wi-[w11 w21 w23]'; 
W2-[w12 w22 w32]'; 
W3-[w13 w23 w33]'; 
W=[W1 W2 W3]'; 
W = W * real(inv(W' * W)^(1/2)); 
wll=W(1,1); 
w12=W(1,2); 
w13=W(1,3); 
w21=W(2,1); 
w22=W(2,2); 
w23=W(2,3); 
w31=W(3,1); 
w32=W(3,2); 
w33=W(3,3); 
for i=1:1:40000 
w1l=sum(((Ww11*data(1,:)+w21*data(2,:)+w3 1 *data(3,:)).3).*data(1,:))/length(data); 
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w21=sum(((w11*data(1,:)+w21*data(2,:)+w3 1 *data(3,:)).3).*data(2,:))/length(data); 
w31=sum(((w11*data(1,:)+w21*data(2,:)+w3 1 *data(3,:)).3).*data(3,:))/length(data); 
w12=sum(((w12*data(1,:)+w22*data(2,:)+w32*data(3,:)).3).*data(1,:))/length(data); 
w22=sum(((w12*data(1,:)+w22*data(2,:)+w32*data(3,:)).3).*data(2,:))/length(data); 
w32=sum(((w12*data(1,:)+w22*data(2,:)+w32*data(3,:)).3).*data(3,:))/length(data); 


w13=sum(((w13*data(1,:) 


w23=sum(((w13*data(1,:)+w23*data(2,:)4 
w33=sum(((w13*data(1,:)+w23*data(2,:)4 


wll-wll-3*wll; 
w12=w12-3*w12; 
w13=w13-3*w13; 
w21=w21-3*w21; 
w22=w22-3 *w22; 
w23=w23-3*w23; 
w3l=w31-3*w31; 
w32=w32-3*w32; 
w33=w33-3*w33; 
Wl=[w11 w21 w23]'; 
W2=[w12 w22 w32]'; 
W3-[w13 w23 w33]'; 
W=[W1 W2 W3]'; 
W = W * real(inv(W' *W)^(1 
Ww112W(1,1); 
w12-W(1,2); 
w13=W(1,3); 
w21=W(2,1); 
w22=W(2,2); 
w23=W(2,3); 
w31=W(3,1); 
w32=W(3,2); 
w33=W(3,3); 
W*W'; 
end 





/2)); 


+w23*data(2,:)+w33*data(3,:)).%3).*data(1,:))/length(data); 


+w33*data(3,:)).%3).*data(2,:))/length(data); 
+w33*data(3,:)).%3).*data(3,:))/length(data); 





W=[w11 w12 w13;w21 w22 w23;w31 w32 w33]; 


y=W'*data; 
figure 
subplot(3,1,1) 
plot(y(1,:)); 
subplot(3,1,2) 
plot(y(2,:)); 
subplot(3,1,3) 
plot(y(3,:)); 
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2. GAUSSIAN MIXTURE MODEL 


Let us consider the data collected from the certain group of people about 
their age and height are represented in two-dimensional vectors. Let the first 
element of the vector be the age and the second element is the height of the 
corresponding person, as described below. P; = [a; hi], Р = [ах №]... Pm = 
[am hm]. Now let us define the problem for classifying the collected group of 
vectors into ‘n’ category called as ‘n’ clusters, such that the centroids of the 
clusters are away from each other and the data belonging to the same cluster 
are nearer to each other as shown in the figure. The dots in the Figure 2-1 
belong to the centroid of the individual clusters. The problem defined is 
called clustering the data, otherwise called as grouping the data. 

Let us model the probability of the particular vector ‘x’ from the 
collected data as the Linear combination of Multi Variate Gaussian density 
function. 

(ie) p(x) = p(ci) * Gaussian density function of ‘x’ with mean vector 
‘т’ and covariance matrix ‘cov,’ + p(c2) * Gaussian density function of ‘x’ 
with mean vector ‘m?’ and covariance matrix ‘cov,’ + ... р(с„) * Gaussian 
density function of ‘x’ with mean vector “m,” and covariance matrix ‘cov,’ 
where p(C,) is the probability of the clustern п. 


(ie) р(х) = р(сі) р(х/сі) + p(c;) р(х/сз) + р(сз) р(х/сз) +. . . p(Cn) 
р(х/с,) 


The model described above is called ав Gaussian Mixture Model. Given 
the data (xi, X2, Xs, X4... Ха), obtaining the mean vectors ту, m», ... and the 
covariance matrices сі, C2, ... and the probability of clusters p(c1), p(c2)... 
such that the probability of the collected data is maximized is the task to be 
performed for modeling. As the collected data are independent to each other, 
the probability of the collected data (р) is represented as the product of p(x1), 


P(X2)..-P(Xm) (ie) Рр-р(хі) р(х›)...р(хь) 
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Figure 2-1. Clustering Example with Four Clusters 


Estimating the best values for mean vectors and covariance matrices such 
that the probability Pp is maximized 15 described below. 

Maximization is obtained by differentiating Pp = Pp=p(x1) p(x2)...P(Xm) 
with respect to the unknown parameters mı, mo, M3, ... My, COV), СОУ», 
COV3...COV, and equate to zero. Solving the above set of obtained equations 
gives the best estimate of the unknown parameters, which are listed below 


m= [(р(су/ху)*(ху) + p(c1/x2) *(x2)+ p(c1/x3) *(x3)+...... p(C1/Xm) *Xm)] 





p(c1/X1)+ p(ci/x2)* p(c1/x3)+...... P(Ci/Xm) 
where p (су/ху) is computed as [p(x1:€1) * р (c))] / PED] 


Note that p(x1C1) is the probability of ‘xı computed using Gaussian 
density function of ‘x’ with mean ‘т’ and covariance ‘cov,’. Also p(x1) = 
p(C1) p(xi/ei) + p (C2) p(x1/C2) + р(сз) p(x1/C3) + ...р(с„) р(ху/с„) 
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The reguirement is to estimate the best value for m1.But for computing 


the best value requires р(с1/х1), which in further requires p(x1/C;) as discussed 
above. p(x1/c1) reguires m1 for computation. [Surprise!]. 


To proceed further, an iteration approach called as Expected — 


Maximization algorithm is used to estimate the best value for m; as 
described below. 


2.1 Expectation-maximization Algorithm 


N — 


. Initialize the mean vectors randomly mi, m», m3,... M, 
. Find the Euclidean distance between the x; vector and the mean vectors 


m, m; .m, are computed as ‘dl’, ‘d2’, ... ‘dn’ respectively. If the 
distance d2 is smaller among all, the x; vector is treated as the vector 
belongs to the 274 cluster. In general if the dk is smallest among all, x, 
vector belongs to the k" cluster. Thus cluster index assigned to the vector 
xı is ‘K’. Similarly cluster index is obtained for the entire data vector. 


. Let the number of data belongs to the first cluster be ‘k’. (ie) Number of 


data vectors having cluster index 1 is ‘k. Then p (cl) = probability of the 
cluster 1 is computed as 


p (C1) = kl 





Total Number of data 
and p(c2) is computed as 
p (C2) = k2 





Total Number of data 
and so on 


. Covariance matrix of the cluster 1 is computed as соу, = E ([X-yx] 


[X-ux')) = E (XX7)-uxux'. For instance, let us consider the vectors x), Хз, 
хт, and хә arranged in the column format, belongs to clusterl and the 
corresponding covariance matrix is computed as follows. 


их = [xi + хз + x; +X9 |4. 


E (XX) = [xi xi xix +хуху + x4 x4! }4 





Thus cov; is computed as E (XX?) - uxux'. 


. Similarly, the covariance matrices cov», covs, соул... COV, are computed. 
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2.1.1 Expectation Stage 


Using the above initialized values, 


[р (c1/X1), p (с/х) p (C1/X3)... p (C1/Xm) 
р (c2/X1), р (C2/X2) p (C2/X3)... p (с/х) 
р (C3/X1), p (C3/X2) p (C3/X3)... p (C3/Xm) 
р (C4/X1), р (C4/X2) р (C4/X3)... р (C4/Xm) 


р (Cn/X1), p (Cn/X2) p (Cn/X3). -P (Cn/Xm)] 
are computed. This is called Expectation stage оЁ Ше E-M algorithm. 


2.1.2 Maximization stage 


Maximization stage of the E-M algorithm belongs to maximizing the 
probability Pp. That is obtained by differentiating the probability Pp with 
respect to unknown parameter and equating them to zero. Solving the set of 
obtained equations gives the best estimate of unknown parameters which are 


listed below. 


[(p(c1/x1)*(X1) + p(ci/x2)*(x5)* р(сі/хз)*(хз)+...... p(C1/Xm) *Xm)] 














m,-— 
p(ci/x1)+ p(c1/x2)+ p(ei/xa)*...... р(с1/х,) 
[(p(ci/xi) *(xi) + p(c1/X2)*(x2)+ р(су/хз)*(хз)+...... P(C1/Xm) *Xm)] 
m= 
p(c1/x1)+ р(су/хә)+ р(су/хз)+...... р(с1/х,) 
[(р(су/х)*(х) + p(ci/x2)*(x5)* p(c1/X3)*(x3)+...... p(C1/Xm) *Xm)] 
з= 
p(ci/xi)* p(ci/x2)* p(ci/xa)t...... p(C1/Xn) 
[(p(c1/x1)*(X1) + p(c1/X2)*(x2)+ p(c1/x3) *(x3)+...... P(C1/Xm) *Xm)] 
Ima4-— 
p(ei/xi)* p(c1/x2)+ p(e1/x3)+...... р(с1/х,) 
[(р(с„/х)®( [х1-их1] [х{-их1]) Tess р(с„/хь)® ( [Xa-Hxn] [Xm-Hxml )] 
СОУ: 








p(c,/x1)* p(Cn/X2)+ р(с,/хз)+...... P(Cr/Xm) 
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The computed values are the current best estimates of means and 
covariance matrices. Using this Expectation stage is executed, followed by 
Maximization stage. Thus Expectation stage and Maximization stage is 
repeated for N iterations and hence best estimate of the means and 
covariance matrices are obtained using E-M algorithm. 

Using the estimated means р(с)), р(с›), p(C3)...p(cn) are obtained and 
hence the Gaussian Mixture Model of the collected is obtained. 


2.2 Example 


Two dimensional vectors are randomly generated and the Gaussian 
Mixture Model of the generated data is obtained using the algorithm 
described above is displayed below. 

About 350 vectors are randomly generated. Among which 100 vectors 
are with mean [0.9 0.8], 100 vectors are with mean [0.7 0.6] and 150 vectors 
are with mean [0.5 0.4]. 

The elements of the individual vector are generated independently and 
hence the estimated covariance matrices are diagonal in nature. 

Estimated values of the GMM model after 10 Iterations are given below 


Mean vector 1 = [0.9124 0.7950] 
Mean vector 2 = [0.7125 0.6091] 
Mean vector 3 = [0.4994 0.3990] 


Covariance matrix 1 
0.0045 -0.0001 
-0.0001 0.0047 


Covariance matrix 2 
0.0048 -0.0003 
-0.0003 0.0048 


Covariance matrix 3 
0.0055 -0.0005 
-0.0005 0.0054 
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2.3 Matlab Program 





gmmmodelgv.m 


%GMM MODEL 

cluster1=clustercol([0.9 0.8],100); 

cluster2-clustercol([0.7 0.6],100); 

cluster3-clustercol([0.5 0.4],150); 

clustertot-[clusterl cluster2 cluster3]; 
Vo[p.q]-kmeans(clustertot',3); This may be used for initializing the means 
a(1,:)=[0.83 0.79]; 

a(2,:)=[0.76 0.74]; 

4(3,:)=[0.6 0.3]; 


%To Estimate q ,covl,cov2,cov3 
[p]=clusterno(clustertot,g); 
[m1,n1 J=find(p==1); 
collect l=clustertot(:,n1); 
if(length(n1)-=0) 
covl=cov(collect1'); 
else 

covl=diag([1 1]); 
end 
[m1,n1 J=find(p==2); 
collect2=clustertot(:,n1); 
cov2=cov(collect2'); 
if(length(n1)-=0) 
cov2=cov(collect2'); 
else 

cov2=diag([1 1]); 
end 
[m1,n1 J=find(p==3); 
collect3=clustertot(:,n1); 
cov3=cov(collect3'); 
if(length(n1)-=0) 
cov3=cov(collect3'); 
else 

cov3=diag([1 1]); 
end 
plot(collect1(1,:),collect1(2,:),'r*") 


hold on 
plot(collect2(1,:),collect2(2,:),'g*') 
hold on 
plot(collect3(1,:),collect3(2,:),'b*") 
hold on 

pause(0.3) 
w1=length(find(p==1))/length(p); 
w2=length(find(p==2))/length(p); 
w3=length(find(p==3))/length(p); 
for j=1:1:10 

pwlx=[]; 

pw2x=[]; 

pw3x=[]; 

for i=1:1:length(clustertot) 


pwlx-[pwlx (w1*pxw(clustertot(:,i),q(1,:),cov1))/... 


(w1*pxw(clustertot(:,i),q(1,:),cov1)+... 





w2*pxw(clustertot(:,1),q(2,:),cov2)+... 
w3*pxw(clustertot(:,1),q(3,:),cov3))]; 


pw2x-[pw2x (w2*pxw(clustertot(:,i),q(2,:),cov2))/... 


(w1*pxw(clustertot(:,i),q(1,:),covl)+... 





w2*pxw(clustertot(:,1),q(2,:),cov2)+... 
w3*pxw(clustertot(:,i),q(3,:),cov3))]; 


pw3x=[pw3x (w3*pxw(clustertot(:,1),q(3,:),cov3))/... 


(w1*pxw(clustertot(:,i),q(1,:),cov1)+... 





w2*pxw(clustertot(:,i),q(2,:),cov2)+.. 
w3*pxw(clustertot(:,i),q(3, соу) 

end 

q(1,1)=sum(pw1x.*clustertot(1,:))/sum(pw1x); 

q(1,2)=sum(pw1x.*clustertot(2,:))/sum(pw 1x); 

q(2,1)=sum(pw2x.*clustertot(1,:))/sum(pw2x); 

q(2,2)=sum(pw2x.*clustertot(2,:))/sum(pw2x); 

q(3,1)=sum(pw3x.*clustertot(1,:))/sum(pw3x); 

q(3,2)=sum(pw3x.*clustertot(2,:))/sum(pw3x); 

cov1(1,1)=sum(pw1x.*[(clustertot(1,:)-q(1,1)).*... 
(clustertot(1,:)-q(1,1))])/sum(pw1x); 

cov 1(1,2)=sum(pw1x.*[(clustertot(1,:)-q(1,1)).*... 
(clustertot(2,:)-q(1,2))])/sum(pw1x); 

cov 1(2,1)=sum(pw1x.*[(clustertot(2,:)-q(1,2)).*... 
(clustertot(1,:)-q(1,1))])/sum(pw1x); 

cov 1(2,2)=sum(pw1x.*[(clustertot(2,:)-q(1,2)).*... 
(clustertot(2,:)-q(1,2))])/sum(pw1x); 
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cov2(1,1)=sum(pw2x.*[(clustertot(1,:)-g(2,1)).*... 
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(clustertot(1,:)-q(2, 1))])/sum(pw2x); 


cov2(1,2)=sum(pw2x.*[(clustertot(1,:)-q(2,1)).*... 


(clustertot(2,:)-q(2,2))])/sum(pw2x); 


cov2(2,1)=sum(pw2x.*[(clustertot(2,:)-q(2,2)).*... 


(clustertot(1,:)-q(2,1))])/sum(pw2x); 


cov2(2,2)=sum(pw2x.*[(clustertot(2,:)-q(2,2)).*... 


(clustertot(2,:)-q(2,2))])/sum(pw2x); 


cov3(1,1)=sum(pw3x.*[( 
(clustertot(1,:)-q(3,1))] 
cov3(1,2)=sum(pw3x.*[( 
(clustertot(2,:)-q(3,2))] 
cov3(2,1)=sum(pw3x.*[( 
(clustertot(1,:)-q(3,1))] 
cov3(2,2)=sum(pw3x.*[( 
(clustertot(2,:)-q(3,2))] 


)/sum(pw3x); 


)/sum(pw3x); 


)/sum(pw3x); 





)/sum(pw3x); 


clustertot(1,:)-q(3,1)).*... 


clustertot(1,:)-q(3,1)).*... 


clustertot(2,:)-q(3,2)).*... 


clustertot(2,:)-q(3,2)).*... 
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[p]=clusterno(clustertot,q); 
w1=length(find(p==1))/length(p); 
w2=length(find(p==2))/length(p); 
w3=length(find(p==3))/length(p); 
[m1,n1 J=find(p==1); 

collect 1=clustertot(:,n1); 

[m1,n1 J=find(p==2); 
collect2=clustertot(:,n1); 

[m1,n1 J=find(p==3); 
collect3=clustertot(:,n1); 

close all 
plot(collect1(1,:),collect1(2,:),'r*") 
hold on 

plot(q(1,1),q(1,2),'co') 
plot(collect2(1,:),collect2(2,:),'g*') 
plot(q(2,1),q(2,2),'mo') 
plot(collect3(1,:),collect3(2,:),'b*") 
plot(q(3,1),q(3,2),'yo') 

pause(0.3) 

totcoll {j}=collect1; 

totcol2 {j}=collect2; 

totcol3 {j}=collect3; 

qcol {j}=q; 

end 
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clustercol.m 
function [res] =clustercol(pos,N) 
xl=rand(1,N)/4-(0.25)/2; 
yl-rand(1,N)/4-(0.25)/2; 
res=[x1+pos(1);y1+pos(2)]; 
clusterno.m 
function [q]=clusterno(clustertot,q) 
d1=(clustertot'-repmat(q(1,:),length(clustertot), 1)).^2; 
dl1=sum(d1"); 
d2=(clustertot'-repmat(q(2,:),length(clustertot), 1)).^2; 
d2=sum(d2"); 
d3=(clustertot'-repmat(q(3,:),length(clustertot), 1)).^2; 
d3=sum(d3'); 
[p,a]J=min([d1;d2;d3]); 
2.4 Program Illustration 
Data to be Clustered Clustering in Iteration 1 
1 T т т т 1 т т 
PT ret 9% {йк | - Cluster 1 
oer * EM * M! il ] ?8 + 
A nimes, t Mean 
Orr to % Bes ^ T uw 1 
t E AAA | Vector 1 
06r * + 41 05 4 
e vt + | 
ost athe t зе 4 25 | » Cluster 2 
% Ku E. ^ | 
ABL M LS je | * Mean 
Р: täh CH | 
оз Ау 4 оз 1 Vector 2 
92 04 os 06 07 98 09 1 1л 02 04 os os 07 98 09 т Cluster 3 
Clustering in Iteration 5 Clustering in Iteration 10 
1 т T T T T T 1 T 
" A » ? Mean 
[ A ш | Vector 3 
öö + Mt na . 
07. p ы 07 4 
os 06 J 
ast 4 1 05 "i 4 
Dar F 04 4 
03- wi 03 4 
02 04 05 06 07 08 09 Li 11 M ол 0.5 0.6 07 08 89 11 
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3. K-MEANS ALGORITHM FOR PATTERN 
RECOGNITION 


If the groups of n-dimensional vectors are given, there is the need to classify 
the vectors into ‘k’ category. The centroid vectors of each group are used to 
represent the identity for ‘k’ category. (i.e) If the new vector is to be 
classified among the ‘k’ category, the distance between the new vector with 
all the centroids are computed. The centroid corresponding to the smallest 
distance is treated as the identified group. The centroids of all the groups are 
obtained using K-means algorithm as described below. 

Consider the set of normalized marks (i.e. the mark ranges from 0 to 1) 
scored by the 100 students in the class for particular subject. Each mark is 
treated as the vector with one-dimensional. There is the need to classify the 
collected marks into 6 groups for allocating 6 grades. This is done using k- 
means algorithm as described below. 


3.1 K-means Algorithm 


Step 1: Initialize the 6 centroids randomly corresponding to 6 grades . 

Step 2: Classify the marks and identify the grade for the individual marks 
using the initialized 6 centroids as described in the section 3. 

Step 3: Find the mean of the marks collected in the particular grade say ‘1’. 
This is treated as the centroid corresponding to the grade 1 used for 
the next iteration. This process is repeated to compute the new sets 
of centroids corresponding to the individual grade. 

Step 4: Go to step 2. 

Step 5: The process involved from Step 2 to Step 4 is considered as the one 
iteration and this process is repeated for finite number of iterations to 
obtain the final centroids corresponding to each grades. 


3.2 Example 
The particular sets of marks (data) are subjected to k-means algorithm Final 


clusters along with clusters obtained at every iteration are displayed below. 
Note that after 6" iteration, clusters are not changed. 
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Figure 2-2. Illustration of K-Means Algorithm 


Figure 2-3. Data along with the final centroids obtained by k-means algorithm 
In the above-mentioned example, the data ranges from 0.0067 to 0.9925 


and the six centroids obtained in the 10" iteration are displayed below. 
0.0743 0.2401 0.4364 0.6115 0.7763 0.9603 (see the Figure 2-3). 


3.3 Matlab Program for the K-means Algorithm 
Applied for the Example given in Section 3.2 


kmeansgv.m 





a=rand(1,100); 

plot(a,zeros(1,100),'*"); 

b=rand(1,6); 

for iter=1:1:10 

hold off 
M-([(a-b(1)).^2;(a-b(2)).^2;(a-b(3)).^2;(a-b(4)).^2;(a-b(5)).^2:(a-b(6)).^2;]; 
[P,CLUSTERNO]=min(M); 

u=[r,'g',b',c','m'y'J; 

figure 
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for 1=1:1:6 

[x,y J=find(CLUSTERN O==i); 
col=[]; 

for k=1:1:length(x) 

col= [col a(x(k),y(k))]; 

end 
plot(col,zeros(1,length(col)),strcat(u(i),'*")) 
hold on 

b(i)=mean(col); 

end 

pause(0.01) 

end 





4. FUZZY K-MEANS ALGORITHM FOR PATTERN 
RECOGNITION 


Consider the problem described in the section 3 for classifying the 
normalized marks into 6 clusters for the assignment of grades. 
In fuzzy k-means technique, fuzzy set theory is used to obtain the optimal 
values of the centroid. 

In this technique the particular vector (In this problem it is the 
normalized mark scored by the student) belongs to all the 6 clusters with 
different membership values. For instant the vector 0.3 belongs to the 
different clusters with different membership values as given below 


Cluster 1 = {0.3 (0.0001)} 
Cluster 2 = {0.3 (0.0448)} 
Cluster 3 = {0.3 (0.0022)} 
Cluster 4 = {0.3 (0.0031)} 
Cluster 5 = {0.3 (0.9492)} 
Cluster 6 = {0.3 (0.0007)} 


The numbers in bold letters are the corresponding membership values. 
(1.е.) 0.3 belongs to the cluster 1 with membership value 0.0001 and belongs 
to the cluster 2 with membership value 0.0448 and so on. Note that sum of 
the membership values is 1. 
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Consider the vector x1 belongs to the cluster 1 whose centroid is c1, then 
(x1-c1)^2 must be minimized. But if the x1 belongs to cluster 1 with the 
membership value m1, then m1*(x1-c1)^2 has to be minimized. Similarly x1 
belongs to the cluster 2 whose centroid 1s c2. Also x1 belongs to the cluster 2 
with membership value m2. Thus the term m2*(x1-c2)^2 also has to be 
minimized. 

Thus the fuzzy k-means problem 15 treated as the optimization technique 
for obtaining the membership values along with the centroids of the 
individual clusters such that 


100 6 

2 ее 
У ХМ,2%(х-сұ) is minimized 
i=1 kel 


М; is membership value of the vector xi belongs in the cluster К. Xi is the 
i^ vector. c, is the К" centroid. 

The algorithm for obtaining the centroids along with membership values 
is displayed below. 


4.1 Fuzzy K-means Algorithm 


Step 1: Intialize the membership values (Mi) randomly. 
Step 2: Compute the centroids of the 6 clusters. 


С1=[Месїог1* (the member ship value of the vector 1 belongs to the 
cluster 1)? (i.e) Mii 2 + Vector2* (the member ship value of the vector 2 
belongs to the cluster 1)? (i.e.) Mi? +... Vectorl00* (the member ship 
value of the vector 100 belongs to the cluster 1 (i.e.) Mi1oo ]/ 


Sum of the squared values of the membership values belonging to the 
cluster 1. 
Similarly the centroids C2, C3, C4, C5 and C6 are obtained. 


Step 3: Update the membership values M using the current values of the 6 
centroids as given below. 


Mi; -1/ [((vector 1-c1)* / (vectorl —c1)?)* + 
((vector 1-c1) / (vectorl c2? + 
((vector 1-c1)? / (vectorl —c3? + 
((vector 1-c1)? / (vectorl c4? + 
((vector 1-c1)* / (vectorl —c5? + 
((vector 1-c1)? / (vectorl —c6?y ] 
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Mas -1/ [((vector 4-c5)? / (vector4 —с1)?)? + 
((vector 4-c5)? / (vector4-c22)? + 
((vector 4-c5)? / (vector4 —c3?)* + 
((vector 4-c5)? / (vector4 —c4°) + 
((vector 4-c5)? / (vector4 —c5°) + 
((vector 4-с5)? / (vector4 —c6°) ] 





Similarly other membership values are computed. 


Step 4: Compute the sum of the sguared difference between the previous 
membership value and the current membership value. If the computed 
value is not less than the threshold value go to step 2 to compute the 
next set of centroids and followed by next set of membership values. If 
the threshold value is less than the threshold value, stop the iteration. 
Thus the centroids are obtained using fuzzy k-means algorithm. Using 
the computed centroids, clustering can be obtained as described in the 
section 3. 


4.2 Example 


The particular sets of marks (data) are subjected to Fuzzy k-means 
algorithm. Final clusters along with clusters obtained at every iteration are 
displayed below. 


Figure 2-4. Illustration of Fuzzy K-means algorithm 
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Figure 2-5. Data and the final centroids obtained using Fuzzy k-means algorithm 


In the above mentioned example data ranges from 0.0099 to 0.9883. 
The final centroids obtained after 10 iteration using fuzzy k-means 
algorithm is displayed below 0.0313 0.2188 0.3927 0.5317 0.6977 
0.8850 (see figure 2-5). Also the graph between the sum of the sguared 
difference between the previous membership value and the current 
membership value (vs.) Iteration number is displayed below. 


50 





Ghange in the membership values 














Iteration 


Figure 2-6. Illustration of change in the membership value in every iteration 
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4.3 Matlab Program for the Fuzzy k-means Algorithm 
Applied for the Example given in the Section 4-2 





fuzzykmeansgv.m 


%Fuzzy k-means algorithm 
a=rand(1,100); 
s=[]; 
fm=[]; 
for i=1:1:100 
r=rand(1,6); 
r-r/sum(r); 
fm=[fm;r]; 
end 
%Computation of centroids using fuzzy membership 


for 1=1:1:10 
for 1=1:1:6 
cen {1} -sum((fm(:,1).^2).*a')/sum(fm(:,1).^2); 
end 
9o Updating fuzzy membership value (fm) 
fml=[]; 
for p=1:1:6 
for g=1:1:100 
temp=0; 
for r=1:1:6 
temp =temp+(((a(q)-cen {p})*2)/((a(q)-(cen {гу ))^2))^2; 
end 
fm1(q,p)=1/temp; 
end 
end 
s=[s sum(sum((fm-fm1).^2))]; 
fm=fm 1; 
b=cell2mat(cen); 
M-([(a-b(1)).^2;(a-b(2)).^2;(a-b(3)).^2;(a-b(4)).^2;(a-b(5)).^2;(a-b(6)).^2;]; 
[P,CLUSTERNO]=min(M); 
u=['r','g','b','c','m','y']; 
figure 
for i=1:1:6 
[x,y J=find(CLUSTERN O==i); 
col=[]; 
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for k=1:1:length(x) 

col= [col a(x(k),y(k))]; 

end 
plot(col,zeros(1,length(col)),strcat(u(i),'*')) 
hold on 

end 

pause(0.01) 

end 

figure 

plot(s) 

xlabel('Iteration') 

ylabel('Change in the membership values 





5. MEAN AND VARIANCE NORMALIZATION 


Suppose in the speech recognition system the two speech signals 
corresponding to the same word are to be compared. Two signals are not 
recorded exactly with the same volume (1.e.) the two signals are not with the 
same energy. Thus before extracting features from the speech signals, there 
is the need to normalize the speech signal in terms of mean and variance so 
that two speech signals are made recorded with the same volume. 

Mean and variance normalization of the signal is obtained as described 
below. Consider the speech signal ‘A(t)’ with mean “па” and variance ‘va 
and speech signal ‘B(t)’ with mean m and variance ‘v’. The requirement is to 
normalize the speech signal B(t) so that the mean and variance are changed 
to the desired mean ‘ma and desired variance ‘vq’ respectively. The 
procedure for the mean and variance normalization is given below. 


5.1 Algorithm 


Step 1: Create the vector completely filled up desired mean ‘mg’. Number of 
elements of the vector is equal to the number of samples in the 
speech signal. 

Step 2: Each and every element of the vector is added with +А or -A. If the 
sample value in the original vector is greater than mean ‘m’, add +A 
to the corresponding sample in the vector created. Otherwise the 
value is added with ‘-A’. The value of the ‘A’ is computed using the 
sample value, mean ‘m’, variance “у” of the signal “A (t)’ and the 
desired variance ‘vq’ . 


Let Д; be the value calculated for the 1" sample of the signal ‘A’ 


(original signal) and is computed as described below. 
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The value А (i) is deviated from “m” with (A (1)-m) when the standard 
deviation of the signal is v? then the deviation of the signal from the 
desired mean “my” with the standard deviation of the signal уд? is given by 


Ai =(va /v) P *(A(1)-m) 
5.2 Example 1 


Consider the speech signal A(t) with mean=0.0062 and уапапсе=0.0057 
which is kept as the reference signal. Consider speech signal 
B(t) corresponding to the same word with mean=0.0061 апа 
variance=0.0094. The speech signal B(t) is subjected to mean and variance 
normalization as described above to obtain B(t) with desired mean and 
variance as shown in the figure below. 


Signal A with mean 0.0062 and variance 0.0057 
T T T T 


06 T t 












































28 1 1 1 1 1 1 1 1 1 
9 0.5 1 1.5 2 25 3 3.5 4 45 X 104 
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0.4 4 
0.2- 4 
0 4 
0.2 4 
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25 І І 1 1 П 1 1 
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05 m 
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Figure 2-7. Illustration of Mean and Variance Normalization of the speech signal 
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Note that the signal at the top and the bottom of the figure 2-7 are almost 
identical. This preprocessing stage is normally used before extracting 


features from the signal for pattern recognition. 


5.3 M-program for Mean and Variance Normalization 





meanvarnormgv.m 


%a is the speech signal to be modified 

%b is the reference signal with desired mean and variance. 
% res is the normalilzed version of signal ‘a’ with desired mean and variance 
m=mean(a); 

v=sqrt(var(a)); 

md=mean(b); 

vd=sgrt(var(b)); 

delta=abs((a-m)*vd/v); 

res=ones(1,length(a))*md; 

[p]-quantiz(a,m); 

p-p5 p=p*2-1; 

res=res+p.*delta; 
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NUMERICAL LINEAR ALGEBRA 
Algorithm Collections 


1. HOTELLING TRANSFORMATION 


Consider the set of n-dimensional vectors collected in which the elements of 
the vector are dependent to each other (1.e.) the covariance matrix computed 
for the collected vectors are not the diagonal matrix. The co-variance matrix 
is computed as follows. 

Let хі, X2,...Xm be the collection of “т”, ‘n-dimensional vectors. The 
co-variance matrix is computed using E[(X-ux) (X-ux) ] = E(XX")- рор, 
where ux is the mean vector computed using the collected vectors (ie) 
um [xi*x?tx37... x,]/m. 

Hotelling Transformation transforms the set of vectors х), X2,...Xm into 
another set of vectors (ie) yi, Y2,...Ym» such that the covariance matrix 
computed for the transformed vectors 1s diagonal in nature. 


Covariance Matrix (CM) aj bi ci di 
b; b; c; d; 
сі €? сз d; 
dı 4 d; dy 
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The element a; in the above matrix gives the information about how the 
first elements of the collected vectors are correlated to each other. Similarly 
сі gives the information about how the first elements are correlated with 
third element. Thus if the matrix is the diagonal matrix, then the elements of 
the vectors collected are made independent to each other. 


1.1 Diagonalization of the Matrix ‘CM’ 


1. Compute the Eigen values and the corresponding Eigen vectors of the 
matrix “СМ? 

Eigen values are computed by solving the determinants | CM-AI| = 0. 
Number of Eigen values obtained is equal to the size of the matrix. For 
instance the size of the above mentioned matrix is 4 and hence number of 
Eigen values = 4. Let it be Ay, А, A3 and Ла 


2. Compute the Eigen vector “X;*corresponding to the Eigen value ‘A,’ 
by solving the equation [CM. - А] X,=0.The size of the Eigen vector ‘X,’ 
is 1X4 for the above example. 


3. Similarly Eigen vectors X2, X; and Х, are computed. 


4. Eigen vectors are arranged rowwise in the matrix to form Transforma- 
tion Matrix (say *TM") 


5. Transformed Vector ‘Y,’ is obtained using the transformation matrix 
as mentioned below. 
Yı= TM*( Хі- u).Similarly Y2,Y3, ...Ym are obtained. If covariance 
matrix is computed for the transformed set of vectors [(1.e.) Yi, Y», Үз... 
Y] the covariance matrix is almost diagonal. 


1.2 Example 


Let us consider the Binary image of size 50x50. Let us collect the positions 
of the Black pixel in the image as the two dimensional vectors. Xposition 
and Yposition are the elements of the vector considered. Hotelling 
transformation is applied as described above to get another set of vectors 
called as transformed positions. The Binary image with all pixels valued “1° 
(i.e.) white is created. Replace the pixel value of the transformed positions in 
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ORIGINAL IMAGE 





TRANSFORMED IMAGE USING HOTELLING TRANSFORMATION 


Dip oe ae en 


—————— 





Figure 3-1. Hotelling transformation for binary image rotation 


the Binary image created with “0° (1.e.) Black. The image obtained is called 
Transformed binary image using Hotelling transformation. Realize that the 
Xposition and Yposition of Black pixels in the Transformed binary image 
are independent to each other. 

Note in the original Binary image, the Xpositions and Ypositions are 
dependent to each other and hence the co-variance matrix 1s not the diagonal 
matrix as mentioned below 


CM- 


102.2500 -80.9500 
-80.9500 102.2500 


But in the transformed Binary image, Xpositions and Ypositions are 
independent to each other and hence the covariance matrix computed is the 
diagonal matrix as given below 


CT = 


21.3000 0.0000 
0.0000 183.2000 


Note that the Diagonal values are Eigen values of the covariance 
matrix C. 
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Note that CM = TM’ * CT * TM 


As the pixel positions of the image are represented as positive integers, 
are transformed vectors are properly biased and rounded off, 


1.3 Matlab Program 





hotellinggv.m 


A=imread('Binimage','bmp'); 

A=double(A); 

[x,y J=find(A==0); 

vector=[x y]; 

Meanvector=sum(vector)/length(vector); 
CM=(vector'*vector)/length(vector); 
CM=CM-Meanvector'*Meanvector; 

%Covariance matrix for the original set of vectors is CM 

[E. V]-eig(C); 

%Transformation Matrix 

TM=E'; 

transvector- TM*(vector' -repmat(Meanvector',1 length(vector))); 
xposition=fix(transvector(1,:)+abs(min(transvector(1,:)))+1); 








yposition=fix(transvector(2,:)+abs(min(transvector(2,:)))+1); 
B=zeros(50,50); 
for i=1:1:length(xposition) 
B(xposition(1),yposition(1))=1; 
end 
%Covariance matrix computed for transformed vectors is computed as CT 
MeanvectorT=sum(transvector')/length(transvector'); 
CT=(transvector*transvector')/length(transvector'); 
CT=CT-MeanvectorT'*MeanvectorT; 
figure 
subplot(2,1,1) 
imshow(A) 
title ORIGINAL IMAGE) 
subplot(2,1,2) 
imshow(uint8(B*255)) 
title TRANSFORMED IMAGE USING HOTELLING TRANSFORMATION") 
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2. EIGEN BASIS 


Group of vectors are described as the vector space. АП the vectors in the 
space are spanned by the particular set of orthonormal basis known as Eigen 
basis as described below. (i.e) Any vector in the space is represented as the 
linear combination of Eigen basis. 


2.1 Example 1 


The two dimensional vector are randomly generated and are plotted as given 
below. The vectors mentioned in the diagram are the Eigen vectors. It is 
computed as follows. 
e Compute the co-variance matrix of the generated two 
dimensional vectors 
e Compute the Eigen values and the corresponding Eigen vectors. 
e Direction of the Eigen vectors computed are parallel to the 
vectors mentioned in the diagram. Note that they are orthonormal 
to each other. 


In the above described example, Eigen vectors are E;= [-0.2459 -0.9693] 
and Е =[-0.9693 0.2459]. 





Figure 3-2. Vector space along with Eigen vectors 
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Also, any vector in the generated vector space is represented as the linear 
combination of El and E2. The co-efficients are obtained as the inner product 
of vector and the corresponding Eigen vector. 


For instance, the vector V;=[20 16] is represented as 
C,E;+C,E, 


C; is obtained as the inner product of E; and V, 
C» is obtained as the inner product of E; and V; 


C; = -20.4269 
С» = -15.4513 


Thus Vi = С,Еү+С»Е» 
= -20.4269 *[-0.2459 -0.9693]+ -15.4513 *[-0.9693 0.2459]. 


Eigen vectors corresponding to the insignificant Eigen values are 
contributing less to the vector representation. In such situation number of 
orthornormal Eigen basis to represent the particular vector is reduced. 

Note that the vectors mentioned in the figure 3-3 is obtained by applying 
KLT transformation to the vector space mentioned in the Fig 3-2. The Eigen 
vectors for this space are E;=[1 0] and E,=[0 1] 








Figure 3-3. Vector space after KLT and the corresponding eigen vectors 


3. Numerical Linear Algebra 93 
3. SINGULAR VALUE DECOMPOSITION (SVD) 


As discussed earlier any real square matrix ‘A’ can be diagonalized using 
eigen matrix ‘E’ (every column vector is the eigen vector of the matrix A") 


A-EDE' 
where ‘D’ is the diagonal matrix filled with Eigen values in the diagonal 


Suppose 1f the matrix “А” is the rectangular matrix of size mxn, then the 
matrix ‘A’ can be represented as the follows 


A-U ^ У". This is called as Singular Value Decomposition. 


АА" is the square matrix with size mxm. Using Eigen decomposition 
AAT=UD,U" 


Similarly A'A of size nxn can be reprersented using Eigen 
decomposition as 


АТА= У р, У" 
IfA-U^VT 
AT =V A Tu! 


Note that ^ and ^! are the same as the matrix is the diagonal matrix. 
AA =П^У VOU 
=U”™'U! [Expected Result] 
Similarly 
KAVA UTU AV. 
= V ^l^ V! [Expected Result] 


The diagonal elements of the matrix ^ is obtained as the positive square 
root of the diagonal elements of the matrix D; or D». Note that the diagonal 
elements of the matrix О, and D» are same except the addition or deletion of 
Zeros. 

Thus the matrix A is represented as the product of the Eigen matrix 
obtained from the matrix AA’, Eigen matrix obtained from the matrix АҒА 
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and the Diagonal matrix ^ obtained as the positive sguare root of the 
elements of the eigen values computed for the matrix AA’ (or) АТА 


(1.е.) The Eigen values are same for both the matrices. 


3.1 Example 


The matrix A is applied to SVD as follows. 
A = 
1 2 3 
4 5 6 
U = 
-0.3863 -0.9224 
-0.9224 0.3863 
V = 


-0.4287 0.8060 0.4082 
-0.5663 0.1124 -0.8165 
-0.7039 -0.5812 0.4082 


A= 


9.5080 0 0 
0 0.7729 0 


U*^*y'- 


1.0000 2.0000 3.0000 
4.0000 5.0000 6.0000 


Eigen values of (A'*A) = Е, 


-0.0000 
0.5973 
90.4027 


Eigen values of (A*A’) = E; 


0.5973 
90.4027 


Note that the diagonal elements of the ^ is obtained as the square root of 
Е or E>. 
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4. PROJECTION MATRIX 
4.1 Projection of the Vector ‘a’ on the Vector ‘b’ 
а 
е 


РЬ b 


Figure 3-4. Projection Illustration 

Consider two vectors a and b. The projection of the vector ‘a’ on ‘b’ is the 
scaled version of the vector ‘b’ (i.e) Pb. The vector e is orthogonal to the 
vector b as shown in the figure 3-3. 

e=a-Pb 

— b! (a-Pb) =0 

= b'a=Pb'b 

=P = (b'a)/ (b'b) 


Therefore the projection of the vector ‘a’ оп ‘b’ is given as bP = [b( b'a)]/ 
[(b'b)]. 


The Projection matrix (In this case, ,it is the single element matrix) is 
given as (bb')/( b! b) = [Рм] 


Thus the projection of the vector ‘a’ оп ‘b’ is given as [Рм] a 
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42 Projection of the Vector on the Plane Described 
by Two Column Vectors of the Matrix *X^ 


Figure 3-5. Projection of the vector on the plane 


The Linear combination of certain number of n-dimensional independent 
vectors forms the vector space. Thus vector space is defined as the space 
spanned by the set of independent vectors. 


Consider the vector ‘a’ projected on the plane described by the vectors 
“ХІ? and ‘x2’. Note that the matrix X =[x1 x2]. The vector ‘p’ can be 
represented as the linear combination of the vector ‘xl’ and ‘x2’ (i.e) 
column vectors of the matrix X as illustrated in the figure 3-4 


p=ul*x1+u2*x2 


p=[X] [ul 
u2 


p=[X][s] 
Also the error vector [a-p] is orthogonal to x1 and x2 as illustrated in the 
figure. 
—(a-Xs)! хі =0 and (a-Xs)! x2 =0 
Jointly represented as [a-Xs]! X = 9 
0 
>a X -(Xs)'X 
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—Xa = X! (Xs) 

—s-(X! X)'X'a 

Therefore the projected vector p-[X][s]-X(X' X)'X'a 

Thus projected matrix Рм is defined as X(X' X)'X? .Also the projected 
vector p is represented as [Рм] a 

In general projection of the vector ‘a’ on the vector space spanned by the 
column vectors of the matrix ‘A’ is given as A(AT A)'ATa 


4.2.1 Example 


Consider the vector space spanned by the column vectors of the matrix 


A={1 3 7 6 
4 6 2 4 
1 7 6 4 
The projection of the vector a= [1 2 3] 7 
on the vector space spanned by the column vectors of the matrix A is 
given using the formula А(А! A)''A" a is given below 
Projection matrix 
Py = А(АТАЈ'А" = 
-0.2813 -0.5625 -0.6563 


-0.1875 1.1250 0.0625 
0.0313 0.1875 1.5313 
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Projected vector is computed as А(А! A)'A™ a 


= |-3.3750 
2.2500 
5.0000 


4.2.2 Example 2 


The column vectors А, В апа С аге displayed below. 
A=[10 2 6 5 9 8 5 0 8 dj 
B=[6 8 9 7 2 4 9 9 4 9| 
C=[43 26 44 29 32 34 35 24 35 32] 


The column vector C is related to the column vectors A and B as the 
linear combination as displayed below C = m*A+ n*B. The reguirement is 
to find the optimal value for the scaling constant m and n. 

If C is in the space spanned by the column vectors of A and B, unique m 
and n values can be computed easily. But if C is not in the space spanned by 
the column vectors of A and B, the constants ‘m’ and ‘n’ are the best 
estimated values such that C'=m*A+n*B is in the space spanned by the 
column vectors А and B. 

Also the mean squared error (ie) E {(C-C’) °] is minimized. This is 
obtained using projection matrix as described below. 

The column vector С” is the projection of the vector “С” on the space 
spanned by the vectors А and B. 

Representing the vectors A and B in the matrix column wise to form the 
matrix ‘P’. 


Thus P=| 10 6 
2 8 
6 9 
5 7 
9 2 
8 4 
5 9 
0 9 
8 4 
4 9 
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The projection matrix is computed as PM=P*inv ((P*P'))*P". Using the 
projection matrix, the С” column vector is obtained as the projection of the 
column vector C as C'-PM*C. 


С- ( 44.3140 
25.6450 
39.9154 
32.0289 
31.4915 
33.4768 
36.9648 
22.2118 
33.4768 
34.0142 


Thus the best estimated values of the variable ‘m’ and ‘n’ are computed 
as 


10 6 |1 (44.3140 

2 8 25.6450 
= 2.9506 
2.4680 


MSE is computed as E [(C-C’) *] = 4.1678. Note that it is not possible to 
find the alternate values for the variables ‘m’ and ‘n’ which gives MSE less 
than the 4.1678. 
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5. ORTHONORMAL VECTORS 


The vectors a, b, c are defined as the orthogonal vectors if their inner product 
matrix is the diagonal matrix. (i.e) dot product of the vector a with itself is 
some constant, whereas the dot product of the vector a with b is 0. 

If the diagonal matrix obtained is the identity matrix, then the vectors are 
called as orthonormal vectors. 


5.1 Gram-Schmidt Orthogonalization Procedure 


Consider the set of independent vectors vl, v2 and v3 which spans the 
particular vector space ‘S’. (i.e) All the vectors in the space ‘S’ are repre- 
sented as the linear combination of the independent vectors (v2, v2, v3). 
They are called basis. 

It is possible to identify the another set of independent vectors 01, 02 and 
03 which spans the space S such that the vectors ol, o2 and 03 аге 
orthonormal to each other. The steps involved in obtaining the orthornormal 
vectors corresponding to the independent vectors are described below. 

Let ‘vl’ and ‘v2’ be the independent column vectors. The vector ‘v1’ is 
treated as the reference vector. (i.e) Normalized ‘v1’ is treated as the one of 
the orthonormal vector ‘ol’= 'vl/|vl||. The orthogonal vector ‘p? is 
computed using the projection of the vector ‘v2’ on the vector ‘v1’. 

From the figure 3-5 p=v2- [(01" v2) / (o1! 01)] ol. The orthonormal 
vector is the normalized form of the vector *p'. 

02 = p/ |р 

Similarly the orthonormal vector “03° corresponding to the vector ‘v3’ is 
obtained as follows. 


v2 


v1 


Figure 3-6. Gram-Schmidt Orthogonalization procedure 
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g=v3- [(01' v3) / (o1! o1)] 01] - [(02' v3) / (о2! o2)] 02] 
03 =q / |lall 

5.2 Example 


Consider the set of independent column vectors v1, v2 and v3 as displayed 
below. 


1 2 4 
vl=|2 у2= |1 | УЗ = |2 
3 4 3 


The corresponding set of orthonormal vectors computed using Gram- 
Schmidt orthogonalization procedure is given below. 


0.2673 0.5203 0.8111 
ol = |0.5345 | 02 =|-0.7804 | о3-| 0.3244 
0.8018 0.3468 -0.4867 


The dot product of the vectors 01, 02 and 03 is displayed in the matrix 
form for illustration. 
1.0000 -0.0000 0.0000 


-0.0000 1.0000 0.0000 
0.0000 0.0000 1.0000 


Note that the matrix obtained is the identity matrix as expected. 
5.3 Need for Orthonormal Basis 


Suppose the vector ‘V’ in the vector space ‘S’ is represented as the linear 
combination of the orthonormal basis computed in the section 5.1 


9 
У=| 6 
14 


The co-efficients are obtained as the inner product of the vector ‘V’ and 
the corresponding orthonormal basis. 
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The co-efficients аге 
сі-У! *о1= 16.8375 
c2=V! *02= 4.8558 
c3=V! *03= 2.4333 


Thus the vector V is represented as 


V = с1*01+с2*02+с3*03 





Computing the co-efficients become easier if the basis are orthonormal to 
each other and hence orthonormal basis is required. 

If adding few more independent vectors increases the number of 
independent vectors spanning the space, the space spanned by these vectors 
also increases. The corresponding orthonormal vectors obtained using Gram- 
Schmidt orthogonalization procedure shows that the orthonormal basis are 
the same as that of the earlier except the inclusion of additional few 
orthogonal basis corresponding to the additional independent vectors. Hence 
the co-efficients of the already existing basis remains the same. 

In the previous example if the vector chosen in the space spanned by the 
vectors v1 and v2. 


6 
| 6 | (зау) 
14 


The vector can be represented as the linear combination of ‘ol’ and 
“027 with the corresponding co-efficients 16.0357 and 3.2950 respectively. 
They are computed using inner product technique. Suppose the same vector 
is represented as the linear combination of ‘ol’ ‘o2’ and “03°, the 
corresponding co-efficients are obtained as 16.0357, 3.2950 and 4.4409е- 
015 respectively. Note that first two co-efficients remains unchanged and 
also note that the value of the third co-efficient is almost insignificant for 
representing the vector. This property is used in compression techniques in 
which the significant co-efficients are stored rather than the vector itself. 
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5.4 M-file for Gram-Schmidt Orthogonalization 
Procedure 





gramgv.m 

%Given the independent vectors as the input,Orthonormal vectors as the 
%output computed using Gram-Schmidt Orthogonalization procedure 
VoInput vectors are arranged rowwise 


function [res]-gramgv(x) 


o{1}=x(1,:)/sqrt(sum(x(1,:).%2)); 
for k=2:1:size(x, 1) 





s=x(k,:); 
for m=1:1:k-1 
s=s - ((o{m}*x(k,:)')/(o{m}*o{m}'))*o{m} ; 
end 
o{k}=s/sqrt(sum(s.%2)); 
end 
res=0; 
6. COMPUTATION OF THE POWERS 
OF THE MATRIX ‘A’ 
Consider the matrix A= | 3 4 |. The matrix A is computed as 


1 2 
described below. 


The matrix “А” can be diagonalized using eigen matrix “E” (every column 
vector is the eigen vector of the matrix *A") as described in the section 3 
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A=EDE' 
Where ‘D’ is the diagonal matrix filled with Eigen values in the diagonal. 


Consider the computation of the matrix А? =A A=EDE' EDE'-EDDE' 
-ED'E' 


In general A? is using ED'?ET 


For the given matrix “A”, eigen matrix “E” is computed as 


E= |0.9315 -0.8421| and the diagonal matrix D is given as follows 
0.3637 0.5393 


D=|4.5616 0 
0 0.4384 


Therefore А!® is computed as ED'”E' = 


5.06471.0e+065 7.90871.0е+065 








1.97721.0e+065 3.08751.0е+065 


Note that if all the eigen values in the Diagonal matrix described above is 
less than 1, the matrix A’ converges to the value zero when n tends to 
infinity. 


7. DETERMINATION ОЕ К!" ELEMENT 
IN THE SEOUENCE 


Let us consider the seguence 0, 1, 2, 3, 6, 11, 20, 37... in which fourth 
element is the summation of three previous numbers. Fifth value is the 
summation of 2, 3 and 4 value of the sequence and so on. 
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The problem is to obtain the value of the 100th element of the above 
mentioned seguence. This is achieved using Eigen decomposition as 
described below. 
Let P, be the value of ће К" position of the sequence. 


Рк = Prat Po +Px-3 


Let us define the matrix Ux = |Р 


Prat 
Pk 
Рк+з 
Uru =| Pre 
Pia 


Ош 4111 | Ux 
100 
010 


a 


Let the matrix A be 111 
101 
010 


Represent the vector Up as the linear combinations of Eigen vectors of 
the matrix A. 

The Eigen values and the corresponding Eigen vectors are computed and 

displayed below. 


Eigen values = [Ал 22241 
= [1.8393 -0.4196 + 0.60631 -0.4196 - 0.60631] 


Eigen vectors arranged in columnwise =[ E1 E2 E3] = 
-0.8503 -0.1412 - 0.37521 -0.1412 + 0.3752i 


-0.4623 -0.3094 + 0.44711 -0.3094 - 0.44711 
-0.2514 0.7374 0.7374 
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Represent the matrix U, as the linear combinations of Eigen vectors. 





2, 
1 = 
-0.8503 -0.1412 - 0.37521 -0.1412 + 0.37521 cl 
-0.4623 -0.3094 + 0.44711 -0.3094 - 0.44711 c2 
-0.2514 0.7374 0.7374 c3 
cl -2.0649 + 0.00001 
c2 |= | -0.3520 + 0.19291 
c3 -0.3520 - 0.19291 
U; = А U, 
ACU: 
Ux = AU, 





О, = cl1*E1+c2*E2+c3*E3 
AU,= c1*A*E1+c2*A*E2+c3*A*E3 
= c1*4,*E1+c2*4,*E2+c3*1,*E3 
Similarly 
U, = AU, = cI*A, *E1+c2*4, *E2+c3*1;*E3 
Therefore Ujo91s computed as 


U100 = 1.0e+026 


5.1220 - 0.00001 
2.7848 - 0.00001 
1.5140 - 0.00001 


Therefore the 100th term in the sequence is given as 1.51401.0e+026 
Let us compute 6" term using the method given above 


3. Numerical Linear Algebra 107 
U,=68.0000 - 0.00001 
37.0000 - 0.00001 
20.0000 - 0.00001 


Therefore 6" term is 20 as expected (see the sequence). 


8. COMPUTATION OF EXPONENTIAL 
OF THE MATRIX 


Let A be the matrix, which can be diagonalizable using Eigen matrix as 
A=EDE'. e^ can be computed using series expansion as given below. 


e^ =I + A +(А2)/2! HA? )/3!+..... 
-[«EDE' HEDE" ЕРЕ!) 2! HEDE” ЕРЕ! EDE")/3!+... 
= ED'E'+EDE' + (ЕР?Е!)/2! +(ED°E")/3! +... 


= E[D? + (D/1! )+(D? /2!) + (D? /3!) +(D* /4!)+...] E' 


= EDE! 


or 


0.5257 -0.8507 -0.6180 0 0.5257 -0.8507 
-0.8507 -0.5257 0 1.6180) 1-0.8507 -0.5257 
Therefore e^ is computed as 


Ee? E! = 
0.5257 -0.8507 атыны 0 0.5257 -0.8507 
-0.8507 -0.5257 0 eee -0.8507 -0.5257 


=| 3.7982 2.0143 
2.0143 1.7839 
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9. SOLVING DIFFERENTIAL EOUATION USING 
EIGEN DECOMPOSITION 


Consider the problem of solving the third order differential equation as 
given below. 


Ou/ ӘР +2 @u/ д2 - ди / дї +u=0 with the initial condition given 
U"(0)23, U'(0)-1, U(0)- -1. 


Representing the above equation in the form of matrix representation. 


V = cle! ! E, 4c2 e?! Е, 4 c3 e?! E, 


Where cl, c2 and c3 are the constants obtained using initial conditions. 
E1, E2 and ЕЗ are the eigen vectors of the matrix A corresponding to the 
eigen values 24,252 respectively. 

The Eigen vectors and the eigen values of the matrix A are displayed 


below. 


Ay =-2.2470 „=0.8019 A3= -0.5550 


E1=[-0.8990 0.4001 -0.1781]" 
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Е2- [-0.4484 -0.5592 -0.6973]' 


E3= [0.2600 -0.4686 0.8443] 


Thus V2| U^" 

U’ | is obtained as follows 

U 
07-сі e??? (-0.8990) + c2. e??? * (20.4484) + c3 е 559% (0.2600) 
U’=cl e??*?* (0.4001) + c2 e99?? t (-0.5592) + c3 e9559?* (0.4686) 


U=cl е22470* (-0.1781) + c2. e?9?? t (0,6973) + c3 e 559?! (0.8443) 


The values of the constants can be solved using the initial conditions 


using the set of equations given below. 


10. 


Let us consider the matrix A=} 2 


U"'(0) =3 =(-0.8990) c1 + (-0.4484) c2. + (0.2600) сз 
U'(0) =1 = (0.4001) c1 + (-0.5592) c2. + (-0.4686) c3 
U(0) = -1 = (-0.1781)c1 + (-0.6973)с2 + (0.8443) сз 


Thus с1= -3.4815, с2= -1.5790, с3- -3.2227 


COMPUTATION OF PSEUDO INVERSE 
OF THE MATRIX 

3» 4 
з 4 5 


Decompose the matrix ‘A’ using Singular Value Decomposition. 
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> 
Il 


[P ]| Q IL R] 
-0.6057 -0.7957 | (8.8839 0 0 
-0.7957 0.6057 0 0.2757 0 


-0.4051 -0.5628 -0.7205 
0.8181 0.1288 -0.5605 
0.4082 -0.8165 0.4082 


> 
Il 


Pseudo inverse matrix А” is given as 
A*= {[R] [Q+][P']} 


[Q] = (1/8.8839 0 


0 1/ 0.2757 
0 0 
= |01126 0 
0 3.6268 
0 0 


ОО" = nt 0 || 
-2.3333 1,8333 
-0.3333 0.3333 


1.6667 -1.6667 


AA'= (1 0 
0 1 
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11. COMPUTATION OF TRANSFORMATION 
MATRICES 


If the vectors in the spacel are spanned by the independent vectors ul, u2, 
u3, ... un. Also if the vectors in the space2 are spanned by the independent 
vectors v1, v2 ...vn . They are called basis of the respective spaces. 

Suppose the vector ‘u’ in the space 1 is mapped to the vector v in the 
space 2. 

(.е) Т (и) = у. The vector ‘v’ can be represented as the linear 
combinations of the independent vectors v1, v2, ...vn. 

Similarly the transformation of the basis vector ul in the spacel can be 
represented as the linear combinations of v1, v2, v3...vn. 


Therefore T (ul) =a;, VI+ а;у2+ а; v3 a14V4 +... аі vn 


Similarly Т (42) =a21 VI+ а››у2+ аз; V3+ agyv4 +... аз, vn 


T (un) =an, у1+ а,у2+ an3 V3+ angv4 +... ann vn 


The co-efficients of the vector vl, v2, v3 .. vn in the above mentioned 
equation are arranged in the vector form and are called co-efficient vector as 
given below. 

The vector [a1; a12 413 214... Ain] is the co-efficient vector associated with 
the vector T(ul) which is in the space 2. Similarly the vector 
[a21 322 аз ал. Arn] is the co-efficient vector associated with the vector 
T(u2). 

Consider the vector ul in the space 1 which can also be represented as 
the linear combinations of the basis vectors ul, u2, u3, ...un as given below. 


ul=1*ul+0*u2+0*u3+....0*un. 





In this case [1 0 0 0 ...0] is the co-efficient vector associated with the 
vector ul in the spacel. 
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As mentioned earlier the any vector ‘u’ in the space 1 and the 
corresponding mapped vector v in the space 2 can also have their associated 
coefficient vectors. Let [p1, p2, p3, ...pn] be the co-efficient vector associated 
with the vector ‘u’ and [ql, q2, q3...qn] be the co-efficient vector associated 
with the vector ‘v’. 

The matrix relating the co-efficient vector of the particular vector in the 
space | and the co-efficient vector of the corresponding mapped vector in 


the space 2 is called transformation matrix and is computed as described 
below. 


(1.е.) 
[pl p2 p3 p4 ...pn]' [TM] = [91 92 93 9... qn]* 


The transformation matrix is obtained using the co-efficient vectors 
computed for the T(u1), T(u2),...T(un) in the space 2 where ul,u2 u3,..un 
are the independent vectors which spans the space 1. 


Thus ТМ = [all a21 a31 a41 a51 ...anl 
al2 a22 a32 a42 а52 ...an2 
al3 a23 a33 a43 a53 ...an3 
al4 a24 a34 a44 a54 ...an4 


aln a2n a3n а4п a5n ...ann 


For example the vector ul with co-efficient vector [1 0 0 0 0 0 0...0] is 


mapped to the vector T(ul) whose co-efficient vector is obtained using TM 
as 


all a21 a31 a41 a51 ...anl 1 а11 
а12 а22 а32 a42 a52 ...an2 0 а12 
а13 а23 а33 a43 a53 ...an3 0 | =| 413 
al4 a24 а34 а44 а54 ...an4 0 al4 
aln a2n a3n a4n asn ...ann 0 aln 


The co-efficient vector obtained is as expected. 
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11.1 Computation of Transformation Matrix 
for the Fourier Transformation 


Consider that the vector space 1 is spanned by the 4 basis as given below. 


vl= v2 v3= v4= 


оо о н 
оо н о 
он оо 
- о о о 


Every vector in the space spanned by the basis vector given above is 
mapped to the vectorin vector space 2 called as Fourier space. Fourier 
transformation of'the basis vectors v1, v2, v3, v4 are given as 

Tvl) =[1 1 1 17 

T(v2) = [1.0000 0- 1.00001 -1.0000 0+ 1.00001 ]' 

T(v3) = [1 -1 1 -1f 

T(v4) = [1.0000 0+ 1.0000i -1.0000 0 - 1.00001 ]" 

The basis vectors of the Fourier space is given as 

fl =(1/2)* [I 1 1 I1] 

f2 = (1/2) * [1.0000 0-1.0000i -1.0000 0+ 1.0000i |! 

B = (1/2)*®[1 -1 1 -1f 

f4 = (1/2)* [1.0000 0+ 1.00001 -1.0000 0-1.000011! 


The transformed vector T(v1) is represented as the linear combination of 
f1, f2, f3, f4. Note that the basis are orthonormal basis 
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The co-efficient vector is obtained as the inner product of the vector 
T(v1) with the corresponding basis. 


[T(v1)] f1" = [2] [Note that the inner product complex numbers ‘a’ 
and ‘b’ is computed as a'b*] 


[Т(у1)] 2° = [0] 
[TVD] fl" = [0] 
[Т(у1)] fl" = [0] 


Therefore the co-efficient vector associated with the vector Т(у1) is 
given as |20001! 


Similarly the co-efficient vector associated with the vector T(v2), T(v3) 
and T(v4) is given as [0 2 0 0]', [0 0 2 0]' and [0 0 0 2]' respectively. 


Thus the transformation matrix for the Fourier transform is given as 


2000 
0200 
0020 
0002 


As the transformation matrix is the diagonal matrix with ‘2’ as the 
diagonal elements, the co-efficient vector for any vector in the Vector space 
1 is ‘coef’ then the co-efficient vector for the corresponding mapped vector 
in the Vector space 2 (Fourier domain) is given be 2*c. 


Thus the Fourier transformation of the vector [2 3 4 1] is given as the 


1 1 1 1 252 10 

1 -1 -1 1 3*2 -24 
(1/2) |1 -1 1 -1 4%2| 2| 2 

1 1 -1 -1 1*2 -2-21 


In the same fashion, DCT, DST the transformation matrix is the diagonal 
matrix and hence the transformation values can be easily obtained using 
simple matrix multiplication as described above. 
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11.2 Basis Co-efficient Transformation 

Consider the 4-dimensional vector space which are spanned by the basis 
ul=[1 0001, u2=[0 1001 and u3=[0 0101 and u4=[0 0 0 1]". Consider 
another set of orthonormal basis which spans the space. 

vl =(1/2)* [1 1 1 1f 

v2 = (1/2) * [1.0000 0- 1.00001 -1.0000 0--1.000011! 

v3 = (1/2)* [1 -1 decl] 

v4 = (1/2)* [1.0000 0+ 1.00001 -1.0000 0- 1.00001 ]' 

Any vector in ће space can be represented as the linear combination of 
ul, 12, 13 and u4. The same vector in the space can be represented as the 
linear combination of v1, v2, v3 and v4. 

Consider the transformation T which transforms the vector v is 


represented as T(v) and is equal to у. 


у is represented as the linear combination of ul u2 u3 and u4.Let the co- 
efficient vector be [pl p2 p3 p4] 


T(v)-v is represented as the linear combination of v1,v2,v3,v4.Let the 
co-efficient vector be [g1 q2 q3 94]. 


The transformation matrix which converts the co-efficient vector 
[pl p2 p3 p4] into the co-efficient vector |41 42 43 q4] is transformation 
matrix for the change of basis. It 1s obtained as described below. 





T([1 0 0 0])=[1 0 0 0] is represented as 0.5%у1--0,5%у2-0,5%у3--0,5%у4 


T([0 1 0 0])=[0 1 0 0] is represented as 
0.5*у1 + (0.5 i )*v2+ (-0.5)*v3+ (-0.5 1)*v4 


T([0 0 1 0])=[0 0 1 0] is represented as 
0.5*v1 + (- 0.5)*v2+ (0.5)*v3+ (- 0.5 ) *v4 


T([0 00 1])=[0 0 0 1] is represented as 
0.5*v1 + (- 0.5 i)*v2+ (- 0.5)*v3+ (0.5 1) *v4 
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There the transformation matrix which converts the co-efficients of the 
particular basis into the co-efficients of the another set of basis belongs to 
the same space is given as 








05 05 05 0.5 
0.5 0.51 -0.5 -0.51 
0.5 -05 0.5 -0.5 
0.5 -0.51 -0.5 0.51 


For instant consider the vector [2 3 4 5] сап be represented using the 
basis ul, 12, u3, u4 with co-efficients [2 3 4 5]. The vector [2 3 4 5] can be 
represented using the basis v1, v2, v3, v4 using the co-efficients computed 
as 


05 05 05 0.5 2 7 
0.5 0.51 -0.5 -0.51 3і|-|-ізі 
0.5 -0.5 0.5 -0.5 4 -1 
0.5 -0.51 -0.5 0.51 5 -I+ 


(i.e.) 7*(v1)H-1-D*v2+(-1)*v3+(-1+)*v4 = [2 3 4 5] 


Thus the transformation matrix which converts the co-efficients of the 
particular vector represented using the basis “17 into the co-efficients of the 
same vector represented using the basis ‘v’. 

If there are only few significant co-efficients obtained when represented 
using the particular set of basis, data compression is achieved. (See chapter 
4). This property is called data compaction property of the transformation. 
Discrete Cosine Transformation (DCT) is having very high data compaction 
property. Hence JPEG (Joint Photographic expert group) is using DCT for 
image compression. Also JPEG2000 standard is using DWT (Discrete 
Wavelet Transformation) for image compression. 


3. Numerical Linear Algebra 117 


11.3 Transformation Matrix for Obtaining Co-efficient 
of Eigen Basis 


Consider the group of two dimensional vectors collected randomly forms the 
vector space. This is consider as the subspace spanned by the independent 
vectors ul = [1 0 ] and u2=[0 1] '. The subspace is spanned by the Eigen 
basis E1=[-0.2459 -0.9693]' and E; =[-0.9693 0.2459]' (see section 2-1). 

The vector [1 0] is represented as the linear combination ul and u2 as [1 
0] = 1*u1+0*u2.The transformed vector of [1 0] is the same vectors itself 
(i.e.)[1 0].The vector [1 0] is represented as the linear combinations of Eigen 
basis given by -0.2459*E1+(-0.9693)*E2 =[ 1 0]. Similarly the vector 
[0 1] is represented as the linear combination of ul and u2 as [0 1] = 
O*ul+1*u2 and its transformed vector [0 1] is represented as the linear 
combination of Е1 and E2 as given by -0.9693*E] +0.2459*E2. 

Thus the transformation matrix which converts the co-efficients of the 
basis ul and u2 into the co-efficients of the Eigen basis Е1 and E2 for the 
particular vector in the space. 


TM = |-0.2459 -0.9693 
-0.9693 -0.2459 


Note that the transformation matrix consists of Eigen basis arranged 
row wise. 


114 Transformation Matrix for Obtaining Co-efficient 
of Wavelet Basis 


Consider the 8-dimensional space spanned by the 8 independent vectors 
ul=[10000000]' u2=[01000000]' u3=[00100000]' 


u4=[00010000]' u5=[00001000]' u6-[00000100] 
u7=[00000010]' u8=[00000001]' 
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Consider any vector in the space can also be represented as the linear 
combination of Wavelet basis as described below. Note that orthonormal 
basis 


wl= (1/sqrt(8)) [11111111] 
w2= (l/sqrt(8)) [1 1 1 1-1-1 -1-1]' 
w3- (1/sgrt(4)) [1 1 -1 -1 00 0 0]' 
w4= (l/sqrt(4)) [O 0 0 0 1 1 -1 -1]* 
w5- (l/sqrt(2)) [1-100000 0]" 
w6- (1/sqrt(2)) [0 0 1 -1000 0]' 
w7= (1/sqrt(2)) [0 0 0 0 1 -10 0]' 
w8- (1/sgrt(2)) [0 00 000 1 -1]' 


As described in the section 11.3 the transformation matrix which 
converts the co-efficient of the basis ul,u2,...u8 into the co-efficient of the 
basis wl, w2, w3,...w8 for the particular vector is computed and is displayed 
below. 













0.3536 0.3536 0.3536 0.3536 0.3536 0.3536 0.3536 0.3536 
0.3536 0.3536 0.3536 0.3536 -0.3536 -0.3536 -0.3536 -0.3536 


0.5000 0.5000 -0.5000 -0.5000 0 0 0 0 
0 0 0 0 0.5000 0.5000 -0.5000 -0.5000 
0.7071 -0.7071 0 0 0 0 0 0 
0 0 0.7071 0.7071 0 0 0 0 
0 0 0 0 0.7071 0.7071 0 0 


0 0 0 0 0 0 0.7071 0.7071 


Note that the transformation matrix consists of the wavelet basis arranged 
row wise. 


12. SYSTEM STABILITY TEST USING 
EIGEN VALUES 


Input signal and output signal of the system are usually related with 
differential equation. The output signal is solved using the eigen 
decomposition as described in the section 9. The general expression of the 
output signal is of the form output (t)=cl е! E, +c2 е?" E; + c3 е! Б, 
where Al, A2 and A3 аге the eigen values. Ej, E2 and ЕЗ are the eigen 
vectors of the matrix described in the section 9. 
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From the above eguation it is observed that Output (t) is bounded only 
when all the eigen values are less than 0 (1.e.) negative values. This is the test 
for stability of the system using eigen values. 


13. POSITIVE DEFINITE MATRIX TEST FOR 
MINIMA LOCATION OF THE FUNCTION 
F (X1, X2, X3, ... XN) 


Consider the function 2 х + 2 x +2 x -2 xi x; -2 хо x3. To find the 
minima of the function f(x), partial differentiate the function with respect to 
х1, x2, x3 and equate to zero and solve for хі, x2, x3 to get (xo, yo, Zo). The 
slope at this point (xo, yo, zo) 1s zero. To confirm that the points so obtained 
are minima, all the partial second derivative values are positive. This is 
tested using matrix technique as described below. 

Represent the equation as the product of three matrices 


[ххх] (2 -1 0 Xi 
-1 2 -l X2 
0 -1 2 X3 


Ifthe matrix | 2 -1 0 | is the positive definite matrix, then the 
-1 2 -I 
0 -1 2 





function f(x) = 2 Xp +2 xj +2 ху^-2 xi x; -2 хо ху is the minima function 
and the minima occurs at the point where the slope is zero (i.e.) first 
derivative is zero. If all the eigen values of the matrix are positive, then the 
matrix is positive definite matrix. Thus sign of all the eigen values of the 
matrix defined above decides whether the point (хо, Yo, Zo) at which the slope 
is zero is minima or not. 


14. WAVELET TRANSFORMATION USING 
MATRIX METHOD 


The 1D signal can be analyzed by using wavelet transformation. This can be 
used to compress the data and to denoise the signal. Wavelet transformation 
of the 1D signal contains the information regarding the low frequency 
content of the signal, which is called approximation co-efficients and the 


120 Chapter 3 
information regarding high frequency content of the signal, which is called 


detail co-efficients. The wavelet transformation of the signal using matrix 
method is described below. 


14.1 Haar Transformation 


Consider the signal with number of samples N. Form the Haar matrix with 
size NXN and with the diagonal matrices filled up with the matrix given 


below. 
1 М, 
1 - 16 


For №8, the matrix obtained is 


^ 0 0 0 0 0 0 
% -2 0 0 0 0 0 0 
00% о 0 0 0 
TMI 00 ”-” 0 0 0 0 
0 0 0 0 ^ ^0 0 
0 0 0 0 4-4 0 0 
000 0 0 0 ^ ^ 
00 0 0 0 0 ^ -% 


Multiply the matrix TM1 with the is the signal data [40 41 42 43 44 4546 
d7 ] to obtain first level decomposition of Haar transformation. 


I level Haar decomposition of the signal 

[2 (40-41) %2 (92-13) % (d4+d5) % (46-47) | 

The samples % (40-41) % (d4+d5) are the average samples. They аге 
called approximated co-efficients of the signal at the first level. The samples 


№ (42-43) % (46-47) are called detail co-efficients of the signal at the first 
level. 


Thus Approximation 1 = [2 (d0+d1) ⁄2(d4+d5)] 


Detail 1- ['4 (42-43) %(46-а7) ] 
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The approximated co-efficients of the signal obtained in the first level 
(Approximation 1) is further decomposed using Haar transformation to 
obtain approximation and detail co-efficients of the second level using the 
transformation matrix as given below. 


^ % 0 0 
ТМ2 = ^ - 0 0 
0 0 ^ ^ 
0 0 ^-^ 


Multiply the matrix [TM2] with the approximated data obtained in the 
first level decomposition to obtain approximation and detail co-efficient of 
the signal at the second level. 


Thus second level approximation is given as 


Approximation 2 = [% ( % (d0+d1) + > (d4+d5) )] 





= %[ d0+d1+d4+d5] 
Detail2= [%(% (d0+d1) - % (d4+d5) )] 


= %[d0+d41- d4 - d5] 


Note that approximation co-efficient in the first level and second level 
are the low frequency information derived from the signal. Similarly the 
detail co-efficient in the first level and the second level are the high 
frequency information derived from the signal. 

The Approximation 2, Detail 1 and Detail 2 completely describes the 
signal. It is possible to reconstruct the signal using Approximation 2, Detail2 
and Detail 2 co-efficients of the signal. 


Reconstruction of the signal is obtained using the inverse Haar matrix 
Formed with the diagonal matrices filled up with the matrix 


m 
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For instance steps involved in reconstructing the signal for N=8 is given 
below. 


е Form the vector with the elements filled up with approximation 2 and 
detail 2 


ГА [ d0+d1+d4+d5] 1⁄4 [d0 +d 1 - d4 - d5]] 


e Form the Inverse transformation matrix 1 


1-00 
1 -1 0 0 
[ITM1] = 0 0 1 1 
0 0 1-1 


e Multiply the vector with the matrix to obtain Approximation 1 and detail 
1 in the jumbled order as Г> (40-41) %2(d2-d3) % (d4+d5) (46-47) | 


e Form the inverse transformation matrix ITM2 


ITM] =Й 1 0 0 0 0 оо 
1-1 000 0 0 0 
0 0 110 0 0 0 
00 oup 0000 
00001100 
00 0 01-1 0 0 
0000001 1 
00000 01-1 


e Multiply the vector obtained in jumbled order as mentioned above with 
the matrix [ITM2] to obtain the original signal [40 41 42 43 44 45 4647] 


1411 Example 


Consider the signal x(n) =a=sin(2*pi*n)+sin(2*pi*100*n) with number 
of samples =512 and sampling frequency Fs =512. The Haar 
transformation is applied to the signal. The approximation and detail co- 
efficients are obtained as described above and is displayed below for 
illustration. Note that approximation co-efficients is the low frequency 
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Figure 3-7. Original signal used for Haar decomposition 
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information derived from the signal and detail co-efficients is the high 


frequency information derived form the signal. 


The signal is reconstructed back using inverse Haar transformation and 


the original signal along with reconstructed signal is displayed below. 
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Figure 3-8. Approximation co-efficients obtained using Haar transformation 
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Figure 3-10. Illustration of Inverse Haar Transformation 
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14.1.2 M-file for haar forward and inverse transformation 





haartrans.m 


i=0:1/511:1; 
a-sin(2*pi*1)*sin(2*pi*100*1); 
plot(a) 
figure 
s=length(i); 
for i=1:1:log2(s) 
temp=createhaarmatrix(length(a))*a'; 
temp=reshape(temp,2,length(temp)/2); 
approx {i}=temp(1,:); 
det {1}=temp(2,:); 
a=approx {i}; 
end 
for k=1:1:(length(approx)-1) 
subplot((length(approx)-1),1,k) 
plot(approx {k}) 
title(strcat('approx',num2str(k))) 
end 
subplot( 
figure 
for k=1:1:(length(approx)-1) 
subplot((length(approx)-1),1,k) 
plot(det[k]) 
title(strcat('detail'Ánum2str(k))) 
end 





createhaarmatrix.m 


function [res]-createhaarmatrix(m) 

temp 1=[ones(1,m/2);(-1)*ones(1,m/2)]; 

temp 1=reshape(temp 1, 1,size(temp1,1)*size(temp1,2)); 
dl=temp1.*(1/2); 

temp2=[ones(1,m/2);zeros(1,m/2)]; 
temp2=reshape(temp2,1,size(temp2,1)*size(temp2,2)); 
d2=(1/2)*temp2(1:1:length(temp2)-1); 
res=diag(d1)+diag(d2,1)+diag(d2,-1); 
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haarinvtrans.m 


app=approx {8}; 

for i=(log2(s)-1):-1:1 
a=[app;det {i} ]; 
a=reshape(a, | ,size(a,1)*size(a,2)); 
app=createinvhaarmatrix(length(a))*a'; 
app=app'; 

end 

figure 

subplot(2,1,1) 

1=0:1/511:1; 

a=sin(2*pi*i)+sin(2*pi* 100*1); 

plot(a) 

title('Original signal"); 

subplot(2,1,2) 

plot(app) 

title'Reconstructed signal); 





createinvhaarmatrix.m 


function [res]=createinvhaarmatrix(m) 

temp 1=[ones(1,m/2);(-1)*ones(1,m/2)]; 

temp 1=reshape(temp 1, 1,size(temp1,1)*size(temp1,2)); 
dl=temp1; 

temp2=[ones(1,m/2);zeros(1,m/2)]; 
temp2=reshape(temp2,1,size(temp2,1)*size(temp2,2)); 
d2=temp2(1:1:length(temp2)-1); 
res=diag(d1)+diag(d2,1)+diag(d2,-1); 
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14.2 Daubechies- 4 Transformation 

Consider the signal consists off N samples. Form the Daubechies 
Transformation Matrix [DTM 1] with diagonal matrices filled up with the 
matrix ‘D’ given below. 


[14N3)/4N2]. [(+V3)/4V2] [G-N3)/4N2] [(1-N3)/ 4\2] 
[((1-V3)/4V2] -[(3-V3)/4V2]  [GN3)/ 4\2] -(1-N3) / 4N2] 


For simplicity matrix D is represented as follows. 
a0 al a2 a3 


b0 b1 b2 b3 


For N=8, the matrix ‘DTM?’ is formed as given below. 


a0 al a2 a3 0 0 0 0 0 0 
bO bl b2 b3 0 0 0 0 0 0 
0 0 a0 al a2 a3 0 0 0 0 
0 0 bO bl b2 b3 0 0 0 0 
0 0 0 0 a0 al a2 a3 0 0 
0 0 0 0 bO bl b2 b3 0 0 
0 0 0 0 0 0 a0 al a2 a3 
0 0 0 0 0 0 bO bl b2 b3 


Note that the matrix ‘D’ are arranged with overlapping in the diagonal of 
the matrix ‘DTM 1’. Compare with the Haar transformation in which the 
matrix are arranged without overlapping. Also note that the size of the 
matrix is of size 8X10 for the signal of size 1x8. So 1x8 sized signal is 
extended to the size 1x10 with 9" and 10" samples are filled up with 1“ and 
274 samples respectively. [Assuming that the signals are repeated cyclically]. 

The steps involved in obtaining the approximation and detail co- 
efficients are as described in the section 14.1 for Haar transformation. 
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Similarly to reconstruct the signal, inverse Daubechies matrix is formed 
with the diagonal matrices filled up with the matrix [ID] as given below. 


[ID] 2| аз b3 al bl 
a4 b4 a2 b2 


For N=8 the second level Inverse Daubechies matrix DITM2 is given as 


a3 b3 al bl 0 0 0 0 0 0 
a4 b4 a2 b2 0 0 0 0 0 0 
0 0 a3 b3 al bl 0 0 0 0 
0 0 a4 b4 a2 b2 0 0 0 0 
0 0 0 0 a3 b3 al Bl 0 0 
0 0 0 0 a4 b4 a2 B2 0 0 
0 0 0 0 0 0 a3 B3 al bl 
0 0 0 0 0 0 a4 B4 a2 b2 


Similar to the DTM 1,size of the matrix is 8x10. Also similar to forward 
transformation, 1x8 sized wavelet co-efficients obtained in the forward 
transformation is extended to 1x10-sized co-efficients with Ist and 2™ co- 
efficients filled up with the last two co-efficients of the wavelet co-efficients. 

The steps involved in reconstructing the signal are same as that of the 
Haar transformation except that in Daubechies matrix the diagonal matrices 
are arranged with overlapping whereas in Haar matrix there is no 
overlapping. 


14.2.1 Example 


Consider the signal x(n) =a=sin (2*pi*n)+sin(2*pi*100*n) with number of 
samples = 128 and sampling frequency Fs = 128. The Daubechies 
transformation is applied to the signal The approximation and detail co- 
efficients are obtained as described above and is displayed below for 
illustration. Note that approximation co-efficients is the low frequency 
information derived from the signal and detail co-efficients is the high 
frequency information derived from the signal. 
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Figure 3-11. Original signal used for Daubechies 4 wavelet decomposition 
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Figure 3-12. Approximation co-efficients at different levels 
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Figure 3-13 Detail Co-efficients at different levels 
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Figure 3-14 Illustration of Inverse Daubechies 4 transformation 


14.2.2 
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M-file for daubechies 4 forward and inverse transformation 
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daub4trans.m 


1=0:1/127:1; 
a-sin(2*pi*1)*sin(2*pi*100*1); 
plot(a) 
figure 
s-length(a); 
for i-1:1:log2(s) -2 
a=[a a(1:2)]; 
i 
temp=createdaubmatrix(length(a)-2)*a'; 
temp=temp(1:1:length(temp)); 
temp=reshape(temp,2,length(temp)/2); 
approx {i}=temp(1,:); 
det {i}=temp(2,:); 
a=approx {i}; 
end 
for k=1:1:(length(approx)-1) 
subplot((length(approx)-1),1,k) 
plot(approx {k}) 
title(strcat('approx',num2str(k))) 
end 
figure 
for k=1:1:(length(approx)-1) 
subplot((length(approx)-1),1,k) 
plot(det[k]) 
title(strcat('detail'Ànum2str(k))) 
end 
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createdaubmatrix.m 


function [res]=createdaubmatrix(m) 


a=[(1+sqrt(3))/(4*sqrt(2)) (3+sqrt(3))/(4*sqrt(2)) (3-sqrt(3))/(4*sqrt(2)) 
(1-sqrt(3))/(4*sqrt(2))]; 
b=[(1-sqrt(3))/(4*sqrt(2)) -(3-sqrt(3))/(4*sqrt(2)) (3+sqrt(3))/(4*sqrt(2)) 


-(1 +sgrt(3)y(4*sgrt(2))]; 

tl=repmat([a(2) b(3)],1,(m/2)); 

t2-repmat([a(3) b(4)],1,m/2); 

t3-repmat([a(4) 0], 1,m/2); 

t4=repmat([b(1) 0],1,m/2); 

res-diag(repmat([a(1) b(2)],1,m/2)) + diag(t1(1:1:m-1),1) 
+diag(t2(1:1:m-2),2)+diag(t3(1:1:m-3),3)+diag(t4(1:1:m-1),-1) ; 

res=[res [zeros(1,m-2) res(1,3) гез(2,3)]' [zeros(1,m-2) res(1,4) res(2,4)]']; 





daubdinvtrans.m 


app-approx {log2(s)-2}; 

for i=(log2(s)-2):-1:1 
a-[app;det {i} ]; 
a=reshape(a, | ,size(a,1)*size(a,2)); 
a=[ a(length(a)-1:1:length(a)) a]; 
app-createdaubinvmatrix(length(a)-2)*a'; 
арр=арр”; 

end 

figure 

subplot(2,1,1) 

1=0:1/127:1; 

a=sin(2*pi*i)+sin(2*pi* 100*1); 

plot(a) 

title('Original signal"); 

subplot(2,1,2) 

plot(app) 

title'Reconstructed signal); 
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createdaubinvmatrix.m 


function [res]=createdaubinvmatrix(m) 


a=[(1+sqrt(3))/(4*sqrt(2)) (3+sqrt(3))/(4*sqrt(2)) (3-sqrt(3))/(4*sqrt(2)) 
(1-sqrt(3))/(4*sqrtQ))]; 
b=[(1-sqrt(3))/(4*sqrt(2)) -(3-sqrt(3))/(A*sqrt(2)) (3+sqrt(3))/(4*sqrt(2)) 


-(1+sqrt(3))/(4*sqrt(2))]; 

tl=repmat([b(3) a(2)],1,(m/2)); 

t2=repmaf([a(1) b(2)],1,m/2); 

t3=repmaf([b(1) 0],1,m/2); 

t4=repmat((a(4) 0],1,m/2); 

res-diag(repmat([a(3) b(4)],1,m/2)) + diag(t1(1:1:m-1),1) 
+diag(t2(1:1:m-2),2)+diag(t3(1:1:m-3),3)+diag(t4(1:1:m-1),-1) ; 

res=[res [zeros(1,m-2) res(1,3) res(2,3)]' [zeros(1,m-2) res(1,4) res(2,4)]']; 
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SELECTED APPLICATIONS 
Algorithm Collections 


1. EAR PATTERN RECOGNITION USING 
EIGEN EAR 


Ear pattern recognition is the process of classifying the unknown ear image 
as one among the finite category. The following is the report on the 
experiment done on ear pattern recognition with the small database. The 
experiment uses twelve ear images collected from four persons. Three 
images are collected from each person. Among them, eight images are used 
to train the classifier. Remaining four images are used to test the classifier. 
The steps involved in Ear pattern recognition using Eigen ears are 
summarized below. 


1.1 Algorithm 


Step 1: Mean and variance of the collected ear images are made almost 
equal using mean and variance normalization technique described in 
the section 3-5. Mean and variance of one of the image from the 
collections is treated as desired mean and desired variance. (See 
figure 4-1) 
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Step 2: 


Step 3: 


Step 4: 


Step 5: 


Step 6: 


Step 7: 


Step 8: 


Step 9: 
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Reshape the matrix into the vector whose elements are collected 
column by column from the matrix. 


Co-variance matrix of the collected vectors is computed. Eigen 
values of the co-variance matrix and the Eigen vectors 
corresponding to the significant Eigen vectors are computed (In this 
example 6 Eigen vectors are computed) 


Eigen vectors are reshaped into the matrix of the original size. They 
are called Eigen ears as given in the figure 4-2 


The Eigen ears are orthogonal to each other. They can be made 
orthonormal to each other by normalizing the vectors. 


For every ear image matrix, feature vectors are obtained as the inner 
product of eigen basis vectors and the reshaped ear image matrix . 


Mean vector of the feature vectors collected from the same person is 
treated as the template assigned to that corresponding person. This is 
repeated for other persons also. Thus one template is assigned to 
every person and they are stored in the database. 


To classify the unknown ear image as one among the four cate- 
gories, template is computed as the inner product of Eigen basis 
vectors (Eigen ears) with reshaped normalized unknown ear image. 
The template thus obtained is compared with the group of templates 
stored in the database using Euclidean distance. 


Template corresponding to the minimum Euclidean distance is 
selected and the person corresponding to that template is declared as 
the identified person. 
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SAMPLE EAR IMAGES BEFORE AND AFTER MEAN.VARIANCE NORMALIZATION 





PERSON 1 AFTER NORMALIZATION PERSON 2 AFTER NORMALIZATION 





PERSON 3 AFTER NORMALIZATION PERSON 4 AFTER NORMALIZATION 


Figure 4-1. Sample Ear Images before and after normalization 


Figure 4-2. Eigen Ears corresponding to the largest Eigen values 
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1.2 M-program for Ear Pattern Recognition 





earpatgv.m 


clear all 
close all 
k=1; 
md=117; 
vd=139; 
for i=1:1:4 
for j=1:1:1 
s=streat('ear' num2str(i),num2str(j),’.bmp'); 
temp=imread(s); 
temp=rgb2 gray(temp); 
temp 1=reshape(temp, 1 ,size(temp,1)*size(temp,2)); 
temp 1=double(temp 1); 
temp 1=meanvarnorm(temp 1,md,vd); 
final {k}=temp 1; 
k=k+1; 
end 
end 
f=cell2mat(final'); 
f=double(f); 
c=cov(f); 
[E,V]=eigs(c); 
for i=1:1:6 
F-reshape(E(:,1),85,60); 
subplot(2,3,1) 
colormap(gray) 
imagesc(F) 
end 


%Feature vector 


for 1=1:1:4 
s=strcat('ear',num2str(1),num2str(1),".bmp'); 
s=imread(s); 
s=rgb2gray(s); 
s=double(s); 
sl=reshape(s,1,size(s,1)*size(s,2)); 


sl=meanvarnorm(s1,md,vd); 
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Е=[]; 
for j=1:1:6 
F=[F sum(s1.*E(:,j))]; 
end 
FEA {i}=F; 
end 
save FEA FEA E md vd 





testinggv.m 


load FEA 
Vo Testing 
[filename, pathname, filterindex] = uigetfile('*.bmp', 'Pick an BMP-file'); 
s-imread(strcat(pathname,filename)); 
s=rgb2gray(s); 
s=double(s); 
sl=reshape(s,1,size(s,1)*size(s,2)); 
sl=meanvarnorm(s1,md,vd); 
Е=[]; 
for j=1:1:6 
F=[F sum(s1.*E(:,j)')]; 
end 
S=[sum((FEA {1 }-F).^2) sum((FEA (2) -F).^2) sum((FEA {3 }-Е).^2) sum((FEA (4j -F).^2)]; 
[P.Q]-min(S); 
switch Q 
case 1 
A-imread('earl 1 bmp"); 
A=rgb2gray(A); 
g=gray(256); 
msgbox('First person','Pattern Recognition','custom',A,g) 
case 2 
A=imread('ear2 1 bmp"); 
A=rgb2gray(A); 
g=gray(256); 
msgbox('Second person','Pattern Recognition','custom', A,g) 
case 3 
A=imread('ear3 1.bmp'); 
A=rgb2gray(A); 
g=gray(256); 
msgbox('Third person','Pattern Recognition','custom',A,g) 
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сазе 4 

A=imread('ear41.bmp'); 

A=rgb2gray(A); 

g=gray(256); 

msgbox( Fourth person','Pattern Recognition','custom',A,g) 
end 





1.3 Program Illustration 


4. Selected Applications 141 


2. EAR IMAGE DATA COMPRESSION USING 
EIGEN BASIS 


Collected ear images as mentioned in the section 1 are represented as the 
group of 85x60 sized matrices. The elements of the matrices are stored with 
the numbers ranging from 0 to 255.8 bits are required to store every number. 
(1.е.) 8*85*60—40800 bits are required to store the single gray image. 
Therefore there is the need to identify the technique for storing the image 
data with reduced number of bits. This technique of representing the data 
with reduced number of bits by removing the redundancy from the data is 
called Image data compression. 

Ear image data compression using eigen basis is the compression 
technique which exploits psycho visual property of the eye. This technique is 
used to compress the similar images belong to the particular group. In the 
experiment described below, the group considered is the ear images 
collected from many persons. The steps involved in Image data compression 
using Eigen basis are described below. 


2.1 Approach 


Step 1: 8x8 sized subblocks of the ear image matrix are collected randomly 
from the ear image database. 


Step 2: Reshape the subblock matrix into the vector of size 1x64. 


Step 3: Thus set of 100 vectors are collected randomly as described in step1 
and step2. 


Step 4: Co-variance matrix is computed for the collected vectors. 


Step 5: Eigen values are computed for the covariance matrix. Eigen vectors 
corresponding to the significant Eigen values are computed. They are 
called Eigen basis, which spans the Ear vector space, which consists of 
all the vectors available as the reshaped sub blocks in the ear images. 
Note that Eigen vectors thus obtained are orthonormal to each other. 


Step 6: Number of Eigen basis obtained as described in the step 4 is less 
than 64 (vector size). This number is called as the dimension of the 
Ear vector space. 


Step 7: To compress the ear image, ear image matrix is divided into 
subblocks of size 8x8. Reshape every sub block into the vector of 
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size 1x64. Represent this vector as the linear combination of Eigen 
basis obtained as described in the step 4. The co-efficients are 
obtained as the product of transformation matrix with the vector. 
Transformation matrix is the matrix filled up with Eigen basis 
arranged in the row wise. Number of Eigen co-efficients obtained is 
egual to the dimension of'the Ear vector space. 


Step 8: Thus every 1x64 sized vector obtained from every sub blocks of the 
ear image is mapped into the vector of size [1 Xb], where ‘b’ is the 
dimension of the ear vector space which is less than 64. Thus every 
sub blocks of the ear image is stored with ‘b’ values which is very 
much less than 64. In this experiment, the value of ‘b’ is 
20<<64. Thus compression with the ratio 3.2:1 is achieved using 
eigen basis technique. 


Step 9: Every sub blocks of the decompressed image are obtained as the 
linear combinations of Eigen basis with the corresponding co- 
efficients associated with that sub block. The compressed and 
decompressed image of the sample ear image is displayed in the 
figure 4-3. 


EAR IMAGE BEFORE COMPRESSION 


CORRESPONDING EAR IMAGE AFTER 
COMPRESSION USING EIGEN BASIS 





Figure 4-3. Ear image compression with the ratio 3.2:1 using Eigen basis 
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2.2 M-program for Ear Image Data Compression 





earcompgv.m 


m-1; 
for i-1:1:4 
for j=1:1:3 

for k-1:1:20 
s=strcat('ear',num2str(1),num2str(7),'.bmp'); 
temp=imread(s); 
temp=rgb2gray(temp); 
x=round(rand*((size(temp, | )-8)))+1; 
y=round(rand*((size(temp,2)-8)))+1; 
t=temp(x:1:(x+7),y:1:(y+7)); 
final {m}=double(reshape(t, 1,64)); 


m=m+1; 


— 





— 


end 
end 
end 
data=cell2mat(final"); 
c=cov(data); 
[E. V]-eig(c); 
% Collecting significant Eigen values 
E-E(:,45:1:64); 
%Transformation matrix 
TM=E'; 
save TMEIG TM 
% Image to be compressed 
A-imread('earl 1.bmp'); 
B=rgb2gray(A); 
C=blkproc(B,[8 8],'compresseig(x)'); 
D=blkproc(C, [size(E,2), 1],'decompresseig(x)); 
figure 
subplot(2,1,1) 
imshow(B) 
title EAR IMAGE BEFORE COMPRESSION?) 
subplot(2,1,2) 
D=D(1:1:size(B, 1), 1:1:size(B,2)); 
imshow(uint8(D)) 
title CORRESPONDING EAR IMAGE AFTER COMRESSION USING EIGEN BASIS") 
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compresseig.m 


function [res]=compresseig(x) 
x=double(reshape(x,1,64)); 
load TMEIG 

res=TM*x'; 
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decompresseig.m 


function [res]=decompresseig(x) 
load TMEIG 
res=reshape(TM'*x,8,8); 
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Figure 4-4. Basis used for Ear Image Compression 
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3. ADAPTIVE NOISE FILTERING USING 
BACKPROPAGATION NEURAL NETWORK 


Consider the signal transmitted through the noisy channel gets corrupted 
with the noise. The corrupted signal is given as the input to the FIR filter as 
given below to filter the noisy part ofthe signal. 

Let x(n) be the noisy corrupted input signal to the system, y(n) be the 
output signal of the system which is the filtered signal and h(n) is the impulse 
response of the system. They are related mathematically as given below. 


y(n) = x(n)*h(n) 


=> 
N-1 

y(n) = > h(k)x(n-k) 
k=0 


‘ 15 the convolution operator. “М” is the order of the filter. For 
instance if the order of the filter is 11, 


y(n) = h(0)x(n) + h(l)x(n-1) + h(2)x(n-2) + h(3)x(n-3) + h(4)x(n-4)* 
-..h(10) x(n-10) 


The h(0), h(1), ... h(10) are the impulse response of the system, otherwise 
called as the filter co-efficients of the system. 

Obtaining the values of the filter co-efficients is the task involved in 
designing the digital filter to do specific operation. To obtain the values 
there is the need to study about the nature of the noise of the channel. 

In practical situations the reference signal is sent through the channel and 
the corresponding corrupted signal obtained as the output of the channel is 
stored. These are used for determining the filter co-efficients of the FIR 
filter. Once filter co-efficients are obtained, they are used to filter the real 
time noisy signal corrupted due to channel transmission at the receiver side. 


Figure 4-5. FIR FILTER 


146 Chapter 4 


FIR Filter block as described above can also be replaced with the Back 
propagation Neural Network block as mentioned in the chapter 1 to perform 
filtering operation as described below. In this experiment the reference 
signal and its corresponding corrupted signal are assumed to be known. 


3.1 Approach 


Step 1: Back propagation neural network (see chapter 1) is constructed with 
11 input Neurons and 1 output neuron and 5 Hidden neurons (say) 


Step 2: In signal processing point of view, input of the neural network is the 
corrupted signal and the output of the neural network is the filtered 


signal. 


Step 3: During the training stage, the elements of the Input vectors are the 
samples collected from the corrupted reference signal. Similarly 
element of the output vector is the corresponding sample collected 
from the reference signal. 





x1 
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x4 
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Figure 4-6. BPNN Filter Structure 
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Table 4-1. Sample Hetero Association table relating input vectors and output vectors 
ofthe BPNN 





х1 х2 х3 х4 х5 хб х7 x8 x9 x10 x11 Y 
X(10 хо) X(8) X(7) X(6) X(5) X(4) XG) X(2) ха) X(0) Y(0) 
X(11) X(10) X(9) X(8) X(7) X(6) X(5) X(4) X(3) X(2) X(1) Ү(1) 
Х(12) X(11) X(10) X(9) X(8) X(7) X(6) X(5) X(4) X(3) X(2) Y(2) 





X(16  X(15 х(4) X03) X02) хап) х(0 хо X) XD Х(6)  Y(16 
ха хау х(5) хал) хаз) X(2) хап) хой хо) X8) XT) YAT) 
X(18 хау) х(6) х(5) х(4) X03 ха? XI) х(0)  X(9)  X(8  Y(18) 
X09  X(18 X7 хп X05) X04 X03 X(2) хп) хао хо) Ү(19) 
хоо X(19 хаз) ха)  X(16 X05) х(4) хаз) X2) XII) хо) Ү(20) 











The Hetero association table relating the sample input vector and the 
corresponding output vector used for training the constructed ANN is given 
below. Note that x is the corrupted reference signal and y is the reference 
signal 


Step 4: Train the Artificial Backpropagation Neural Network as described 
in the chapter 1. Store the Weights and Bias. 


Step 5: Now the BPNN filter is ready to filter the original corrupted signal. 


3.2 M-file for Noise Filtering Using ANN 





noisefiltuanngv.m 


%Noise filtering using ANN 
A=wavread('C:\Program Files\NetMeeting\Testsnd'); 
B=A(:,2); 

B=B(6001:1:7000,1); 

B=B'; 

B=(B+1)/2; 
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%Assumed channel characteristics which is disturbing the signal. 


H=[0.7071 -0.7071]; 
C=conv(H,B); 


%Collections 


INPUTVECTOR= ]; 
OUTPUTVECTOR=| ]; 


for i=1:1:5000 
x=round(rand*(100-12))+1; 
INPUTVECTOR=[INPUTVECTOR;C(x:1:x+10)]; 
OUTPUTVECTOR=[OUTPUTVECTOR;B((x+11))]; 
end 


W 1=rand(5,11); 

В1=гапа(5,1); 

W2=rand(1,5); 

B2=rand(1,1); 

nntwarn off 

[W1,B1,W2,B2]-trainbpx 
(W1,B1,'logsig', W2,B2,'logsig', NPUTVECTOR',OUTPUTVECTOR',[250 100000 0.02 
0.01]; 

save WEIGHTSNOISEFILT W1 ВІ W2 B2 

D=blkproc(C,[1 11,0 5],'dofilter(x)'); 

D=(D-0.5)*2; 

figure 

subplot(3, 1,1) 

plot(B) 

title( Original signal") 

subplot(3,1,2) 

plot(C) 

title('Corrupted signal") 

subplot(3,1,3) 

plot(D) 

title('Retrieved signal") 
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dofilter.m 


function [res]=dofilter(x) 
load WEIGHTSNOISEFILT 
res=simuff(x',W1,B1,'logsig',W2,B2,'logsig’); 





3.3 Program Illustration 


Note that low frequency information is lost and that is not retrieved in the 
filtered signal. This is due to the fact that first 100 samples are used to train 
the ANN. If the training sequence is increased, low frequency information 
can also be preserved in the filtered signal. Also note that MATLAB Neural 
Network toolbox is used to train the network. 


Original signal 
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Figure 4-7 Noise Filtering using BPNN 
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4. BINARY IMAGE ROTATION USING 
TRANSFORMATION MATRIX 


СҰ 
[707 0.7071 | 9707, 0.7071 | 


450 





m 


Figure 4-8. Vector basis for Binary Image Rotation 


Consider the vector space spanned by the basis [1 0]' and [0 1]! . 
Transformation of the vector [1 0] is the vector [0.7071 0.7071] that 1s in 
the same space. Transformed vector is represented as the linear 
combination of the basis [1 0]' and [0 1]' with vector co-efficient 0.7071 
and 0.7071 respectively. Similarly the transformation of the vector [0 111 is 
the vector [-0.7071 0.7071] that is in the same space. 

The transformed vector [-0.7071 0.7071] is represented as the linear 
combination of the basis [1 0]' and [0 1]' with co-efficient 0.7071 апа 
0.7071 respectively. 

The transformation described above rotates the vector counter clockwise 
by 45° as mentioned in the figure 4-7. 

Thus the transformation matrix used to rotate the vector counter 
clockwise by 45? is computed using the method described in the chapter 3 is 


given below. 
0.7071 0.7071 
0.7071 -0.7071 
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4.1 Algorithm 

Step 1: Read the Binary image. 

Step 2: Find the position vectors of the black pixels of the Binary image. 


Step 3: Transform the position vectors using the transformation matrix to 
obtain the new sets of positions. 


Step 4: Create the blank image completely filled up with white pixels. Fill 
the selected pixels of the blank image described by the positions 


obtained in step 3 with black to obtain the rotated image. 


ORIGINAL IMAGE 


IMAGE AFTER 45 degree ROTATION 





Figure 4-9. Binary Image Rotation 
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4.2 M-program for Binary Image Rotation 
with 45 Degree Anticlockwise Direction 





binimrotationgv.m 


A=imread(selectfig8',"'bmp'); 

TM=[0.7071 0.7071;-0.7071 0.7071]; 

[X Y]=find(A==0); 

VECT=[X Y]; 

TRANSVECT=TM*VECT'; 

TRANSVECT=round(TRANS VECT-min(min(TRANSVECT))+1); 

m=max(max(TRANSVECT)); 

B=ones (m, m); 

for k=1:1:length(TRANSVECT) 
B(TRANSVECT(1,k), TRANS VECT(2,k))=0; 

end 

figure 

subplot(2,1,1) 

imshow(A) 

title ORIGINAL IMAGE) 

subplot(2,1,2) 

imshow(B) 

title IMAGE AFTER 45 degree ROTATION") 





5. CLUSTERING TEXTURE IMAGES USING 
K-MEANS ALGORITHM 


To retrieve the similar images from the database, sample image looking like 
the same is used as the key. The low level features extracted from the key 
image is used as the index for searching. So the data base has to be 
maintained based on the low level feature vectors derived from the images. 
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Consider 23 texture images are stored in the database. They are indexed 
using the low level features collected from the image itself, so that 23 texture 
images are stored under 4 categories. K-means algorithm is used to classify 
the collected texture into 4 categories so that each image in the database is 
associated with the number indicating the category. (1.е.) 1 indicates first 
category, 2 indicates second category and so on. This is called texture 
grouping and is done as described below. 


5.1 Approach 


Step 1: Every texture image in the database is divided into sub blocks of 
size 8x8. The variance of every sub blocks are calculated. They are 
arranged in the vector form of size 1x64. This is called Low-level 
feature vector extracted from the image itself. Note that elements of 
the low level feature vector can also be any other statistical 
measurements measured from the image. 


Step 2: The collected Low-level feature vectors are subjected to k-means 
algorithm to classify the images into four categories as described in 
the chapter 2 


Step 3: Thus the index number is identified for every image in the database. 
Also the centroids of the individual category are also stored. The 
centroid is the vector which is of the same size as that of the feature 
vector. [Refer chapter 2] 


Step 4: To retrieve the images from the indexed database that looks like the 
image which is given as the key image is described below. 


e Extract the low level feature vector from the key image as described in 
the step 1. 

e Compute the Euclidean distance between the computed feature vector 
and all the centroids corresponding to the individual category. The 
category corresponding to smallest Euclidean distance is selected. 

e All the images in that category is retrieved using the index number 
associated with every images 
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Figure 4-10. Texture images - Cluster 1 


Figure 4-11. Texture images - Cluster 2 


Figure 4-12. Texture images - Cluster 3 


Figure 4-13. Texture images - Cluster 4 
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5.2 M - program for Texture Images Clustering 





patclasgv.m 


for 1=1:1:23 
S=strcat(num2str(1),'.bmp'); 

temp 1=rgb2gray(imread(S)); 

A ti)=temp1(1:1:40,1:1:48); 
temp2=blkproc(A {i},[8 8],'var2(x)'); 
FEA {i}=reshape(temp2, 1,size(temp2,1)*size(temp2,2)); 
end 

FEA VECT=cell2mat(FEA’); 
P=kmeans(FEA VECT,4); 

for 1=1:1:4 

figure 

[x,y]=find(P==1); 

for j=1:1:length(x) 
subplot(1,length(x),j) 

imshow(A {x(j)}) 

end 


end 





var2.m 


function [res ]=var2(x) 
res=var(reshape(x,1,64)); 
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6. SEARCH ENGINE USING INTERACTIVE 
GENETIC ALGORITHM 


Consider the small Image database stored with Low-level feature vectors 
associated with every image. Low-level feature vector associated with every 
image is the 1D vector whose elements are the statistical measurements like 
variance, kurtosis, and higher order moments etc., computed from the sub 
blocks of the image. Low-level feature vector associated with every image in 
the database acts as the index for searching the particular image from the 
database. Interactive Genetic algorithm based searching procedure is 
described below. 


6.1 Procedure 


1. The front panel of the search engine consists of finite number of images 
displayed. (Say 8 images) as shown in the figure given below. 

2. The user is asked to assign the rank for the individual images. The value 
ranges from 0 to100 based on the user's interest. (1.e.) User assigns higher 


Figure 4-14. Sample front panel 
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rank if user likes that image and expecting the similar images in the next 
iteration. 

3. 4 Images are selected from the 8 images using Roulette wheel selection 
(Refer chapter 1) as described below. 


e Let the rank values assigned for the eight images are represented as the 
vector as given below. 


[1090101010 90 10 90] 


Normalized vector is computed as 


[0.0313 0.2813 0.0313 0.0313 0.03133 0.2813 0.0313 0.2813] 


Simulate Roulette wheel 


Cumulative summation of the above-normalized vector is displayed below. 
[0.0313 0.3125 0.3438 0.3750 0.4063 0.6875 0.7188 1] 


Random number is generated as ‘r’. If r is in between ‘a’ and ‘b’ in the 
cumulative summation vector song corresponding to the value ‘b’ is selected 
as mentioned below. 


Generated random number = 0.9501(say) 
Therefore select 8" image. 


Generated number =0.2311(say) 
Therefore select 2" image. 


Similarly. 
0).6068----------------------------------- 6" image 
0.0185----------------------------------- 1° image 


4. Thus the selected images are 8,2,6,5 


5. Low level features corresponding to the image number 8,2 6 and 5 are 
computed. The vectors thus obtained are randomly chosen and is 
subjected to cross over and mutation operation as described below. 


Low level feature (LLF) vector before and after crossover 


Vector 1 =[abcdef] 
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Vector2=[ABCDEF] (say) 


Note that the number of elements shown in the vector is considered as 6. 
[Actually the size of the vector is around 200] 


Low level feature (LLF) vector after crossover 


Vector l'- [abcd EF] 
Vector 2° = [A B C D e f] 


6. Thus 8 Low level feature vectors for the second generation are obtained 
as mentioned below. 
LLF1, LLF2, LLF3, LLF4, LLF5, LLF6, LLF7 and LLF8 (say) 


7. Compute the Euclidean distance between the feature vector LLF1 with 
all the low level feature vectors in the database. Retrieve the image 
corresponding to the LLF vector in the database whose Euclidean 
distance is smallest. Repeat the procedure for the feature vectors 
LLF2, LLF3, ... LLF8 


8. Thus the image selected for the next iteration is listed below. 
9, 4, 16, 3, 13, 18, 45, 43 (say). Note that 9 in the selected list indicates 
9" image in the database. 


9. Repeat the step from 1 to 8 until the user is satisfied with the images 
obtained in the latest iteration. 


6.2 Example 


Consider the Image database consists of 94 Natural sceneries. Every image 
is truncated to the standard size 200x200. The feature vector for the 
particular image is computed as described below. 

The image is divided into sub blocks of size 50x50 with overlapping size 
of 8x8.The first feature namely variance of the histogram of the hue content, 
are computed for the all the overlapping sub blocks of the image. This is 
treated as first part of the feature vector corresponding to the image. 
Similarly other features are computed for all the overlapping sub blocks of 
the image to obtain other parts of the feature vector. All the parts belonging 
to the individual features are combined to obtain the feature vector 
corresponding to that particular image. 

Thus feature vectors are computed for all the images in the data base. 
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6.2.1 List of low level features computed in every sub blocks 
of the image 


1. Variance of the histogram of the hue content of the particular subblock is 
computed. 


2. Measurement of uniformity 


p(Zi is the probability of gray value zi in the particular sub block of the 
image, The uniformity of sub block of the image is measured using the 


formula given below. : 
2 
U= > p (2) 


i-l 


3. Measurement of Average Entropy 


Average entropy is computed as described below. 


ға 
g-- рі) ов; plz) 
1-0 


4. Measurement of Relative smoothness 


1 


+2202) 


where, o(z) is the standard deviation of the gray values 2 in the 
particular sub block of the image. 


5, 6, 7, 8. Variance of the approximation Wavelet co-efficient 


The image is subjected to discrete wavelet transformation (DWT) to 
obtain approximation co-efficient (app), horizontal detail co-efficient (hor), 
vertical detail co-efficient (ver) and diagonal (both horizontal and vertical) 
detail co-efficient(dia) The variance of the approximation co-efficient ‘app’ 
is computed and considered as 5" feature. Similarly the variance of the detail 
co-efficient ‘hor’,’ver’ and ‘dia’ are considered as the 6", 7 and 8" features 
respectively. 
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9. Skewness 


The Skew ness of the gray values is computed for every sub block of 
the image as mentioned below. 


1-1 4 
(2) => (2; - т) р(д) 


i=l 
10. Kurtosis 


The Kurtosis of the gray values is computed for every sub block of the 
image as mentioned below. 


2-1 4 
и (т) = i 7m) plz) 


i=l 


6.3 M-program for Interactive Genetic Algorithm 





igagv.m 


r=round(rand(1,8)*94); 

figure 

for k=1:1:8 

subplot(4,2,k) 

title('Assign the rank') 

s-strcat(num2str(r(k)),' pg"); 
A=imread(strcat('E:\Naturepix\',s)); 

imshow(A) 

title(strcat('Image',num2str(k))); 

end 

for iter=1:1:10 

s1=input('Enter the rank for image 1(Between 0 to 100)'); 
s2=inpuf( Enter the rank for image 2(Between 0 to 100)'); 
s3=inpuf( Enter the rank for image 3(Between 0 to 100)'); 
s4=input('Enter the rank for image 4(Between 0 to 100)'); 
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s5-input('Enter the rank for image 5(Between 0 to 100)'); 
s6=input('Enter the rank for image 6(Between 0 to 100); 
s7-input('Enter the rank for image 7(Between 0 to 100); 
s8-input('Enter the rank for image 8(Between 0 to 100); 
vect=[s1 s2 s3 54 s5 56 s7 s8]/100; 
normvect=vect/sum(vect); 
c=cumsum(normvect); 
u=[ J; 
for 1=1:1:4 
rl=rand; 
с2=с-г1; 
[x,y]=find(c2>0); 
if(y(1)—1) 
ul=1 
else 
ul=y(1)-1 
end 
u=[u ul]; 
end 


load NATUREFEATURE 

for 1=1:1:4 

feature (1) =FEATURE(r(u(1)),:); 

end 

%Cross over to obtain 8 feature vectors 

k=1; 

for i=1:1:4 

pl=round(rand*159)+1 

r=round((rand*3)+1); 

dl=feature fr); 

r=round((rand*3)+1); 

d2=feature fr) 

nextfeature(kj-[d1(1:1:p1) 42(р1+1:1:160)]; 

k=k+1; 

nextfeature(k)=[d1(1:1:p1) 42(р1+1:1:160)]; 

k=k+1; 

end 

г=[]; 

for 1=1:1:8 
n=nextfeature (1); 
d=(FEATURE-repmat(n,94,1)).%2; 
s=sum(d') 
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{m1 ,n1]=min(s); 
r=[r 01]; 


end 


for k=1:1:8 

subplot(4,2,k) 

title('Assign the rank") 
s=streat(num2str(r(k)),' іре); 
A=imread(strcat('ENaturepixY,s)); 
imshow(A) 
title(strcat('Image',num2str(k))); 
end 

end 





igaimagegv.m 


for 1=1:1:94 
s=strcat(num2str(1),". jpg"); 
A=imread(strcat('E:aturepixY,s)); 
В {i}=naturefeature(A); 

end 

FEA TURE=cell2mat(B'); 

save NATUREFEATURE FEATURE 





naturefeature.m 


function [res ]=naturefeature(x) 
x=x(100:1:size(x, 1), 1:1:size(x,2)-100,:,:); 
x=imresize(x,[200 200]); 

y=rgb2hsv(x); 

h=y(:,:,1); 

res=blkproc(h,[32 321,18 8],'huevar(x)'); 
res |=reshape(res, 1 ,size(res, 1) *size(res,2)); 
y=rgb2gray(x); 
res=blkproc(y,[32 321,18 8],'unif(x)'); 
res2=reshape(res, 1,size(res,1)*size(res,2)); 
res=blkproc(y,[32 321,18 8],'entropy(x)'); 
res3=reshape(res, 1,size(res,1)*size(res,2)); 
res=blkproc(y,[32 32],[8 8],'relativesmooth(x)'); 
res4=reshape(res, 1,size(res,1)*size(res,2)); 
res=blkproc(y,[32 32],[8 8],'wavelet(x)'); 
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res5=reshape(res,1,size(res,1)*size(res,2)); 
res=blkproc(y,[32 321,18 8],'skew(x)'); 
res6=reshape(res, | ,size(res,1)*size(res,2)); 
res=blkproc(y,[32 321,18 8],'kurt(x)'); 
res7=reshape(res, | ,size(res,1)*size(res,2)); 
res-[resl res2 res3 res4 геѕ5 res6 res? |; 





huevar.m 


function [res]=huevar(x) 
x=hist(reshape(x,1,size(x,1)*size(x,2),500); 


res-var(x); 





unif.m 


function [res]-unif(x) 
z=reshape(x, 1,size(x,1)*size(x,2)); 
h=hist(double(z),256); 
h-h/sum(h); 

res-sum(h.^2); 





entropy.m 


function [res]-entropy(x) 
z=reshape(x, 1,size(x,1)*size(x,2)); 
h-hist(double(z),256); 
h-h/sum(h); 

res-(-1)*sum(h.*1og2 (h+eps)); 





relativesmooth.m 


function [res]-relativesmooth(x) 
z=reshape(x, 1,size(x,1)*size(x,2)); 
res=1-(1/(1+var(z))); 
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wavelet.m 


function [res]-wavelet(x) 

x=double(x); 

[c,1]=wavedec2(x,1,'db6'); 
a=appcoef2(c,1,'db6',1); 
vl=var(reshape(a,1,size(a,1)*size(a,2))); 
dl=detcoef2('h',c,1,1); 
v2=var(reshape(d1,1,size(d1,1)*size(d1,2))); 
d2=detcoef2('V',c,1,1); 
v3=var(reshape(d2,1,size(d2,1)*size(d2,2))); 
d3=detcoef2('d',c,1,1); 
v4=var(reshape(d3,1,size(d3,1)*size(d3,2))); 
res=[v1 v2 v3 v4]; 





skew.m 


function [res]=skew(x) 
x=double(x); 

z=reshape(x, 1,size(x,1)*size(x,2)); 
h=hist(double(z),256); 
h-h/sum(h); 

m-mean(z); 
res-sum(((z-m).^3).*h(z-1)); 





kurt.m 


function [res]=kurt(x) 
x=double(x); 

z=reshape(x, 1,size(x,1)*size(x,2)); 
h=hist(double(z),256); 
h=h/sum(h); 

m=mean(z); 
res=sum(((z-m).”4).*h(z+1)); 
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6.4 Program Illustration 





166 Chapter 4 


7. SPEECH SIGNAL SEPARATION AND DENOISING 
USING INDEPENDENT COMPONENT ANALYSIS 


Consider the two speakers talking in the meeting simultaneously with no 
background noise. Two microphones are used for recording. One kept very 
close to the first person. Another microphone is kept near to the second 
person. The recorded signals of both the microphones can be treated as the 
linear combinations of independent sources. [Two speech signals] Hence 
ICA algorithm can be used to separate the two independent signals [Refer 
chapter 2]. 

Consider the second situation in which single speaker is talking with the 
background noise. Two microphones are used to record the signal. One 
microphone kept near to the speaker. Another microphone kept near to the 
assumed noise source. The recorded signals of both the microphones can be 
treated as the linear combinations of independent sources [Speech signal + 
Noise ].Similar to the above ICA algorithm can be used to separate the two 
independent signals. Hence the ICA algorithm can be used to separate the 
two independent signals and hence denoising is achieved. 


7.1 Experiment 1 


Two speech signals х1(4) and x2(t) are linearly mixed to get two mixed 
signals y1(t) and y2(t) as given below. 


y1(t)=0.7*x1(t) +0.3*x2(t) 
y2(t)=0.3*x 1(t)+0.7*x2(t) 


The mixed signals are subjected to ICA algorithm. The independent 
signals x1(t) and x2(t) are obtained as shown below. Note that FASTICA 
toolbox is used to run the ICA algorithm. 


Source signal 

















Figure 4-15. Original speech signals 
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Figure 4-16. Mixed signals 
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Figure 4-17 Separated speech signals using ICA 


7.2 Experiment 2 


Speech signal x(t) and noise n(t) are linearly mixed to get two mixed signals 
у1(0 and y2(t) as given below. 


у) = 0.7*x(t) + 0.3*n(t) 
y2(t)= 0.3*x(t)+ 0.7*n(t) 


The mixed signals are subjected to ICA algorithm to obtain the speech 
signal and noise signal separately. 








speech signal 
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Figure 4-20. Retrieved speech and noise signals using ICA 


4. Selected Applications 169 


7.3 M-program for Denoising 





noiseicagv.m 


a=wavread('C:\MATLAB7\work\audio\revsam\ram 1.wav'); 
b=rand(length(a), 1); 

m=min(length(a),length(b)); 

a=a(1:1:m); 

b=b(1:1:m); 

mixed1=0.7*a+0.3*b; 

mixed2=0.3*a+0.7*b; 


mixed=[mixed1'; mixed2']; 





[signal]-fastica(mixed); 
figure 

subplot(2,1,1) 
plot(mixed1) 
title(Mixed signal1') 
subplot(2,1,2) 
plot(mixed2) 
title(Mixed signal2') 
figure 

subplot(2,1,1) 

plot(a) 

title('speech signal") 
subplot(2,1,2) 

plot(b) 

title(noise signal") 
figure 

subplot(2,1,1) 
plot(signal(1,:)) 
title(Retrieved speech signal obtained using ICA") 
subplot(2,1,2) 
plot(signal(2,:)) 
title(Retrived noise signal obtained using ICA") 


Note that the signal b(n) in the above program is another speech signal in 
case of speech signal separation. The rests of the program remains same. 
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8. DETECTING PHOTOREALISTIC IMAGES 
USING INDEPENDENT COMPONENT 
ANALYSIS BASIS 


The digital images captured using the camera are called as photographic 
images and the images which are not captured using the camera are called as 


photorealistic images. ICA basis are used in classifying the produced image 
into photo graphic or photo realistic images as described below. 


Figure 4-21. Sample photo realistic images 


Figure 4-22. Sample photographic images 


4. Selected Applications 171 
8.1 Approach 


e The 10 photographic images and the 10 photo realistic images are 
collected (see figure 4-3 and figure 4-4). The sub blocks of the image 
sized 16x16 are randomly collected from both the categories. 500 sub 
blocks are collected from photographic images and 500 subblocks are 
collected from the photorealistic images. 


e The subblock thus obtained is reshaped into the size 1x256. Auto 
Regressive (AR) co-efficients of size 1x10 are obtained from the 
reshaped subblock. This is repeated for all the collected sub blocks. The 
set of AR vectors collected from photographic images and the photo 
realistic images are subjected to ICA analysis and 10 ICA basis each of 
size 1x10 are obtained. [Refer chapter 2] 


е Every АК co-efficient vector obtained from the particular subblock of the 
images is represented as the linear combination of 10 corresponding ICA 
Basis. The co-efficient of the ICA basis are obtained using the inner 
product of ICA basis with the corresponding AR co-efficients. This is 
called feature e vector of that particular subblock of the image. 


e Thus 500 feature vectors are collected from photographic images and the 
500 feature vectors are collected from the photorealistic images. The 
centroid of the feature vectors collected from the photographic images is 
computed as Cl. Similarly the centroid of the feature vectors collected 
from the photo realistic images are computed as C2. 


8.1.1 To classify the new image into one among the photographic 
or photorealistic image 


The image is divided into sub blocks. Feature vectors are extracted form 
every sub blocks using ICA basis as described above. The Euclidean 
distance between the feature vector obtained from the particular subblock 
and the centroids “СІ? and ‘C2’ are computed as “417 and ‘d2’ respectively. 
Assign the number 1 to that particular subblock if ‘d1’ is lowest Otherwise 2 
is assigned to that subblock. This is repeated for all the sub blocks of the 
image to be classified. 

Count the number of ‘1’s and the number of ‘0’s assigned to the 
subblocks of the image to be classified. If the number of 175 is greater than 
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0’s, decide that the image is photographic image. Otherwise the image is 


classified as photo realistic image. 


8.2 M-program for Detecting Photo Realistic Images 
Using ICA Basis 





photorealisticdetgv.m 


k=1; 
for 1=1:1:20 
s=streat('C:\Documents and Settings\user!\Desktop\Photorealistic\',num2str(i),' jpg’); 
for j=1:1:50 
temp=imread(s); 
temp=rgb2gray(temp); 
temp=temp(101:1:200,101:1:200); 
rl=round(rand*(size(temp, 1)-8))+1; 
r2=round(rand*(size(temp,2)-8))+1; 
temp=reshape(temp(r1:1:r1+7,12:1:12+7),1,64);; 
temp=double(temp); 
A{k}=Ipc(temp, 16); 
k=k+1; 
end 





end 


temp=cell2mat(A"); 
temp=temp(:,2:1:17); 
I=fastica(temp); 
save I I 


k=1; 
for i=1:1:20 

s=streat('C:\Documents and Settings\user!\Desktop\Photorealistic\',num2str(i),' jpg); 
for j=1:1:100 

temp=imread(s); 

temp=rgb2gray(temp); 

temp=temp(101:1:200,101:1:200); 

rl=round(rand*(size(temp, 1)-8))+1; 

r2=round(rand*(size(temp,2)-8))+1; 
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temp=reshape(temp(r1:1:11+7,12:1:12+7),1,64); 
temp=double(temp); 
temp=lpe(temp,16); 
temp=temp(2:1:17); 
a=[ J; 
for n=1:1:size(I,1) 
a= [а sum(temp.*I(n,:)) ]; 
end 
feaextl {k}=a; 
k=k+1; 
end 

end 


FEA-cell2mat(feaext1"); 
Cl=mean(FEA(1:1:1000,:)); 
save СІ СІ 


C2=mean(FEA(1001:1:2000,:)); 
save C2 C2 





phototest.m 


for i=1:1:20 
s=streat('C:\Documents and Settings\user!\Desktop\Photorealistic\',num2str(i),' jpg"); 
A=imread(s); 
B=rgb2gray(A); 
B=B(101:1:300,101:1:300); 
imshow(B) 
C=blkproc(B,[8 8],'assigncenno(x)'); 
if(length(find(C==1))>length(find(C==0))) 
disp(The image is photorealistic image") 
else 
disp(The image is photographic image") 
end 
end 
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assigncenno.m 


function [res]=assigncenno(x) 
x=double(x); 
load Cl 
load C2 
load I 
x=reshape(x, 1,64); 
p=lpe(x,16);; 
temp=p(2:1:17); 
a-[]; 
for n=1:1:size(1,1) 

а= [a sum(temp.*I(n,:)) ]; 
end 
dl=sum((C1-a)."2); 
d2-sum((C2-a).^2); 
if(d1<d2) 

res=0; 
else 

res=1; 


end 





8.3 Program Illustration 
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9. BINARY IMAGE WATERMARKING USING 
WAVELET DOMAIN OF THE AUDIO SIGNAL 


The process of hiding the information like text, binary image, audio etc. into 
another signal source like image, audio etc. is called watermarking. The 
approach involved in watermarking the binary image signal in the wavelet 
domain of the audio signal is described in the example given below. 


9.1 Example 
Step 1: Consider the audio signal sampled with the sampling rate of 11025 Hz. 


Step 2: The audio signal is divided into frames with 1024 samples. Decompose 
every frame of the speech signal using ‘db4’ wavelet transformation 
[See wavelet transformation in chapter 3]. Hide the first bit of the 
binary data collected from the binary image into the third level detail 
co-efficients obtained using wavelet decomposition of the first audio 
frame. This is done by changing the values of the third level detail co- 
efficients such that mean of the third level detail co-efficients is 
modified to ‘m+k’ if the binary value is ‘1’ or to ‘m-k’ if the binary 
value is ‘0’.Note that variance remains same.[See Mean and variance 
normalization in Chapter 2]. The variable “m” is the mean of the fifth 
level detail co-efficients. The value for the constant ‘k’ is chosen as 
1/5" of the energy of the third level detail co-efficients. 


Step 3: This is repeated for the other frames of the speech signal for hiding 
the entire bits obtained from the binary image. 


Step 4: Thus Binary image is hidden in the wavelet domain of the audio 
signal. 


Step 5: The hidden binary data is retrieved by comparing the mean of the 
corresponding detail co-efficients computed before and after hiding. 
If the mean of the third level detail co-efficient of the particular 
frame corresponding to the original signal is greater than the mean 
of the third level detail co-efficients of the respective frame 
corresponding to the watermarked signal, the bit stored is 
considered as ‘1’,otherwise ‘0’. 


Step 6: Note that the duration of the audio signal is chosen such that all the 
binary data collected from the binary image are hidden in the audio 
signal. 
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In this example length of the audio signal used for hiding binary image 
data is 11488 samples. It is repeated 100 times so that the length of the audio 
signal becomes 1148800 samples. The size of the Binary image used is 
45X45.Binary image is stored at the rate of | bit in 256 samples of the audio 
signal. 


9.2 M-file for Binary Image Watermarking 
in the Wavelet Domain of the Audio Signal 





audiowatermarkdb4gv.m 


A= wavread ('C:\WINDOWS\Media\Microsoft Office 2000\glass.wav'); 

A-repmat(A, 100,1); 

А-А”; 

len=fix(length(A)/256); 

for k=1:1:len-1 
temp=A((k-1)*256+1:1:(k-1)*256+256); 
col {k}=daub4water(temp); 

end 

I=imread((BINIMAGE. bmp’); 

I-1(1:1:45,1:1:45); 

I-reshape(L1,45*45); 


for k-1:1:length(1) 
switch I(k) 
case 0 
temp-col {k} (2) {3}; 
e-sum(temp.^2)/5; 
m=mean(temp); 
v=var(temp); 
col {k} {2} {3}=meanvarnorm(temp,(m-e),v); 
case 1 
temp=col {k} {2} {3}; 
e-sum(temp.^2)/5; 
m=mean(temp); 
v=var(temp); 
col{k} {2} {3}=meanvarnorm(temp,(m+e),v); 
end 
end 
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s=256; 

signal=[]; 

for k-1:1:length(1) 
approx=col {k} {1}; 

det=col {k} {2}; 

app=approx {log2(s)-2}; 

for i=(log2(s)-2):-1:1 
a-[app;det {i} ]; 
a=reshape(a, | ,size(a,1)*size(a,2)); 
a=[ a(length(a)-1:1:length(a)) a]; 
app=createdaubinvmatrix(length(a)-2)*a'; 
app=app'; 

end 

signal=[signal арр]; 

end 

figure 

subplot(3,1,1) 

plot(signal(1:1:5520)); 

title( Original signal") 

subplot(3,1,2) 

plot(A(1:1:5520),'r'); 

title('Watermarked signal"); 

subplot(3,1,3) 

plot(signal(1:1:5520)-A(1:1:5520)) 

title('Difference signal") 


save signal signal 





daub4water.m 


function [res]=daub4water(x) 

а=х; 

s=length(a); 

for i-1:1:log2(s) -2 
a-[a a(1:2)]; 

temp-createdaubmatrix(length(a)-2)*a'; 

temp-temp(1:1:length(temp)); 
temp=reshape(temp,2,length(temp)/2); 
approx {i}=temp(1,:); 
det {1}=temp(2,:); 


a=approx {i}; 
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епа 
res {1 }=арргох; 
res {2}=det; 





gethiddenimage.m 


% size of the Binary hidden image is 45x45 
A= wavread ('C:\WINDOWS\Media\Microsoft Office 2000\glass.wav'); 
A=repmat(A,100,1); 
А=А', 
len=fix(length(A)/256); 
for k=1:1:len-1 
temp=A((k-1)*256+1:1:(k-1)*256+256); 
coll {k}=daub4water(temp); 
end 


load signal 

A=signal; 

len=fix(length(A)/256); 

for k=1:1:len 
temp=A((k-1)*256+1:1:(k-1)*256+256); 
col2 {k}=daub4water(temp); 

end 


temp=[]; 
for k-1:1:45*45 
ml-mean(coll {k} (2) {3}); 
m2=mean(col2 {k} {2} {3}); 
if(m1>m2) 
temp=[temp 0]; 
else 
temp=[temp 1]; 
end 
end 
temp 1=reshape(temp,45,45); 
figure 
imshow(temp 1); 
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meanvarnorm.m 


function [res]-meanvarnorm(a,md,vd) 
m=mean(a); 

v=sqrt(var(a)); 

vd=sqrt(vd); 

delta=abs((a-m)*vd/v); 
res=ones(1,length(a))*md; 
[p]-quantiz(a,m); 

р=р/; 

р=р*2-1; 

res=res+p.*delta; 





createdaubmatrix.m 


function [res]-createdaubmatrix(m) 


a=[(1+sqrt(3))/(4*sqrt(2)) (3+sqrt(3))/(4*sqrt(2)) (3-sqrt(3))/(4*sqrt(2)) ... 
(1-sqrt(3))/(4*sqrt(2))]; 


b=[(1-sqrt(3))/(4*sqrt(2)) -(3-sqrt(3))/(4*sqrt(2)) (3+sgrt(3)) /(4*sqrt(2)) ... 
-(1+sqrt(3))/(4*sqrt(2))]; 


tl=repmaf([a(2) b(3)],1,(m/2)); 

t2-repmat([a(3) b(4)],1,m/2); 

t3-repmat([a(4) 0],1,m/2); 

t4=repmat([b(1) 0], 1,m/2); 

res-diag(repmat([a(1) b(2)],1,m/2)) + diag(t1(1:1:m-1),1)+... 
diag(t2(1:1:m-2),2)+diag(t3(1:1:m-3),3)+diag(t4(1:1:m-1),-1) ; 


res=[res [zeros(1,m-2) res(1,3) res(2,3)]' [zeros(1,m-2) res(1,4) res(2,4)]']; 
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createdaubinvmatrix.m 


function [res]=createdaubinvmatrix(m) 

a=[(1+sqrt(3))/(4*sqrt(2)) (3+sqrt(3))/(4*sqrt(2)) (3-sqrt(3))/(4*sqrt(2)) (1-sqrt(3))/(A*sqrt(2))]; 
b=[(1-sqrt(3))/(4*sqrt(2)) -(3-sqrt(3))/(4*sqrt(2)) (3+sqrt(3))/(4*sqrt(2)) -(1+sqrt(3))/(4*sqrt(2))]; 
tl=repmat([b(3) a(2)],1,(m/2)); 

t2-repmat([a(1) b(2)],1,m/2); 

t3=repmat([b(1) 0],1,m/2); 

t4=repmat([a(4) 0],1,m/2); 

res-diag(repmat([a(3) b(4)],1,m/2)) + diag(t1(1:1:m-1),1)+ ... 
diag(t2(1:1:m-2),2)+diag(t3(1:1:m-3),3)+diag(t4(1:1:m-1),-1) ; 


res=[res [zeros(1,m-2) res(1,3) res(2,3)]' [zeros(1,m-2) res(1,4) res(2,4)]']; 





9.3 Program Illustration 


Original signal 








ü 1000 2000 3000 4000 5000 6000 
Watermarked signal 





ü 1000 2000 3000 4000 5000 6000 
Difference signal 











-1 1 1 1 1 1 
ü 1000 2000 3000 4000 5000 6000 


Figure 4-23. Audio signal before and after watermarking 
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Figure 4-24. Original and watermarked signal in zoomed level 


Hidden image Retrieved hidden image 
used for hiding ұрт the watermarked image 


Figure 4-25. Hidden and Retrieved Binary image 
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